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THE STRESS DISTRIBUTION IN THE HEAD OF A 
THIN-WALLED PRESSURE VESSEL 


By H. B. FERGUSSON, J. KUDAR (G. A. Harvey & Co. (London), Ltd.) 
and R. B. HARVEY (Royal Naval College, Greenwich) 


[Received 20 February 1951. Resubmitted 8 November 1951] 


SUMMARY 


The stress distribution in the head of «. thin-walled vessel formed of zones from 
various surfaces of revolution is considered. The equations of the membrane theory 
are written in terms of the meridional are length, and are then shown to have an 
approximate solution valid for moderate thicknesses of the shell, provided that 
R, cot 6 < R,, where R,, R, are the principal radii of curvature and ¢ the inclina- 
tion of the normal to the axis of the shell. The method is applied to a toro-spherical 
vessel head with a large bunged opening in the centre. The maximum stresses in 
this head occur in a hyperboloidal zone, and are found to be considerably less than 
would arise in a comparable toroidal section. 


1, Introduction and summary of previous work 

THE present paper had its origin in an order made to G. A. Harvey & Co., 
Ltd., by the Shell Petroleum Co., Ltd., for four vessels of the form shown 
in Figs. 1 and 2. In order to test the vessels with air to 25 lb./sq. in. it was 
considered necessary to design them for a maximum stress not exceeding 
8 tons/sq. in. at any point when the internal pressure was 45 |b./sq. in. 

A first study by Fergusson and Kudar led to the form shown in Fig. 1 for 
the transition between the 16-ft. and 10-ft. cylinders, in which the small 
cylinder is joined by a hyperboloid to the spherical portion of the head. 
For the ‘discontinuity components’ of the stresses, trigonometrical solutions 
were used. Harvey later replaced these by the present exponential expres- 
sions, thereby simplifying the calculations. 

Under the condition that the bending moments make a negligible con- 
tribution to the resistance to pressure, the membrane stresses in axially 
symmetric vessels are given by Timoshenko (1, p. 457), together with 
differential equations relating these to the meridional and circumferential 
bending moments and to the radial shearing force (normal to the shell). 
These quantities are all determined by the thickness of the shell and its 
principal radii of curvature. Methods of solution of these differential 
equations are given in power series by Watts and Burrows (2), Salet and 
Barthélemy (3), and de Leiris and Barthélemy (4); for very thin shells, 
approximate methods have also been found by Geckeler (5) and Coates (6). 


(Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
5092.21 B 
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For a shell of thickness ¢, whose principal radii of curvature are ?,, in the 
meridian plane, and R,, the power series solution is valid (7) provided that 
t/R, > 1/50 and tR,/ R 


slowly, and for the hyperboloid to be considered these solutions do not apply, 


1/200. For thin shells the series converge very 
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arrangement of the 
pressure vessel. 


Fic. 1 Meridional section of the pressure vessel. 


In Geckeler’s solution, on the other hand, the lower derivatives are 
neglected in the differential equations, and this is justifiable only if 


Ry t/R? <[12(1—v2)] = 3:3, 


where v = Poisson’s ratio. Since at the junction of the hyperboloid and 


the small cylinder R,t/R? = 1-1, this solution also does not appear to be 
useful in the present case. 

Fortunately a comparative study of Geckeler’s and Coates’s methods 
reveals an interesting aspect of the underlying mathematical problem. It 
will be shown in section 2 that Geckeler’s approximation in neglecting d/d¢ 
compared with d?/d¢? is unnecessary,} and that a more exact solution is 


t See section 2 for definition of ¢. 
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obtainable without complicating the treatment. This follows Coates in 
using derivatives with respect to the are length s of a meridional section, 
but instead of aiming at a simple differential equation for the normal 
deflexion, we return to the basic equations which determine the meridional 
rotation and normal shearing force. A suitable form of solution for these 
equations is found for the hyperbolic section in section 3. 

Since all stresses depend on the radii of curvature of the surface, account 
must be taken of the discontinuity of R, at the joins of the various portions 
of the surface. The material unity of the shell, however, prevents discon- 
tinuous expansion on either side of the boundary, and results in radial 
shearing stresses and additional bending moments at the junctions. These 
udditional stresses are referred to as the discontinuity stresses in sections 
1,5. The continuity conditions at the boundaries which are derived directly 


from Hooke’s law have then to be applied to the continuous components. 


2. Transformation of equations 

\xes are taken with the zy-plane as a meridional plane of the cylinder, 
the y-axis being the axis of the pressure vessel. The tangent to the section 
f the shell in this plane makes an angle ¢ with the x-axis. The following 


notation will be used for other quantities: 


p = internal pressure in the vessel. 
t = thickness of shell. 
R, = meridional radius of curvature of shell, 


radius of curvature of section in xy-plane. 
R, = radius of curvature perpendicular to xy-plane, 


' 
xr Ccosec o. 


Q shearing force per unit length normal to surface. 


U OR, 
V = angle of rotation of a surface element of shell. 
Biv Young’s modulus, Poisson’s ratio. 


0}, 0, = meridional and circumferential (hoop) normal stresses. 


s are length of section in xy-plane. 
u meridional deflexion of the shell. 
w normal deflexion of the shell 


Assuming that bending is negligible, the membrane theory gives correct 
values for o, and Oo, namely 
Oo »R,| 2t. Oo o,(2 R,/R ). l 
eat 2 1\ aaa 
Expressions for the normal stresses and differential equations for the- 


bending moments and shearing stresses in the shell are given by Timoshenko 
1, p. 457). 
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Again, assuming uniform thickness t, we have 


L(U)+—U = Ev. é 
(U) R, J (2) 
RO SW tea ete 
(Y) R, D (3) 
where D = Et 12(1 v?), 
R, d? Ltd fk. cot2d 
L = — — + - = =} + ot - . 
Ride R, nlg)+R ‘cot | R, “) 





In the approximate methods of solution, Cini (5) neglects all but the 
first term in L, whereas the simplification used by Coates (6) is in terms 
of derivatives taken along the meridian, so that the are length s satisfies 
ds = R,dd. We shall first consider the effect of this change of variable on 
the operator L for surfaces of the forms appearing in the vessel required. 

We may write L in the form 


R, d/ d (R, cot?d 
L “(r, ~)\+|R, ‘\+5 + Pact oot 
R, ds ( a) | 1 ds ei tl R, 
1? R,dk, dk, Rk, aR _cot! ~ 
ee 4 —3_ —3 — ta ot 
* ds* iz ds ' ds R, ds — cot | 
Now R, sind = « and differentiating with sibs to s we have 
dk, . dd : dx 
7. sin¢d+ R,cos¢ ae 
or a sin +5 2cos¢ = cos ¢. (5) 
‘ ty 
Substituting for dR,/ds, L will reduce to 
1? cot! a) 
L t ees 6) 
R, 7 - co $o- R, (6 
Since ds = R, dd, a measure of the ratio of the coefficient of the first term 


to that of the second, and of the second term to the third, is R, cot ¢d/R,. 
For the spherical portion of the shell, R, = R, = R = constant, 
and with the dimensions in Fig. 1, ¢ > 28°, cotd < 1-88. Further, 
R,t/R? = t/R = 1/227, and is negligible compared with ,/{12(1—v*)} = 3°. 
We can therefore adopt Geckeler’s approximation in this region. 
Taking L = R, d?/ds? leads to the differential equation 


p+ st J i b (7) 


where X may be U or V; this is in any case a known approximate equation 
when X is the normal deflexion w. In (7), v?/R3 has been neglected in com- 


parison with v?/¢?. 
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STRESS 





In the toroidal shell, R, = constant, and denoting by a the distance of 


the corresponding centre from the axis of the toroid, we have 


R, a cosec ¢- R,. (8) 
Thus Reotd/R, ?, /a |, and for this section also we may use 
L = R, d?/ds®; again R,t/R} = 1/7, which is much less than ,/{12(1—v?)}. 


For an ellipsoidal head of semi-axes (a,b), writing k = a/b, we have 
R,/R, 1+ (k?—1)sin?¢ cos*¢ +k? sind. (9) 
Thus R, t/ R? k?t/R,, and R, cotd/R, > cot d/k?, and the approximation 
2 1 1 l 2 
L R,d* ds? can be used again for the region where cot db k?., 


In the hyperbolic section to be considered, taking 


— y" 
= = i, (10) 
gG- b2 
ind writing k alb, v | (k?+-1)sin?d 1} ‘1 the radii of curvature are 


R, = akv, R, akv’. 
Numerically, at the junction with the 10-ft. cylinder, we have ¢ 


oe 
4 


R, 6in., anda 60 in.; thus k? 10. 

In this section the ratio R, cot ¢/R, increases with decreasing ¢, from 
at the cylinder to about 1 near the junction with the spherical zone. 
[he criterion R,t/ R} for the validity of Geckeler’s approximation, on the 
ther hand, varies from about 1/100 at the spherical junction to nearly 
| next to the cylinder. Hence, for one reason or other, the simplified 
perator R, d*/ds* can be used for L in (2), (3) throughout the hyperboloid 
region. 

Lastly, the terms containing v/R, in the differential equations (2), (3) 

ive to be discussed. Eliminating U’ or V from these equations we obtain 


the two fourth-order equations 


12.1—v?)  v ly 
LLM L{ U —0, 11 
| Ry R, | _ 
j 12(1—v?) vr ( v : 
LLV 4 — i, J 0. 12 
| { R? i)| _ 


In these, we may neglect v?/ R? in comparison with 12(1—v*)/t? provided 
that R,/t > v/,/{12(1—v?)} 1/10, a condition which is necessarily satis- 
fied, the least value of R,/t being about 9. Further, Z(1/R,) can be shown 
to reduce to 

3(k?+-1) 
R, R, 


cot?d 


[5 cos*d(v? + 1)—v? sin’ | — RR.” 
1 fte 


ind thus varies between } at the cylinder-hyperboloid junction and about 
02 at the junction with the sphere. These will be neglected in comparison 
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with 12(1—v?)/#2 = 20. We may therefore use for both U and V solutions 
of the differential equation 


9 
LL(X) 4.12! 


l—y?) , 
—_— X= 0. (13) 
t* 

Since the coefficient of X is constant, this equation may be reduced to 
two second-order equations (5, p. 8), (1, p. 458) 


LX+cX = 0, (14) 
where c? 12(1—v?)/é?. 
These two equations may be written as 
d2X were 
— + 2iA2X 0, (15) 
ds? ~~ 


with 22 y[8(1—v?)]/ gt. 

For a constant value of A, the equation (12) has the four independent 
solutions xX o(t1+irs (15a) 
with all four combinations of sign. In the general case of a variable J, i.e. 
variable R,, a solution by the W.K.B. method (7) is suitable when 

(dA/ds) < X. 
For the hyperbola considered, 
A? const/ R, 
so that I 


: const dR, rv’ il dk. l l 
2AA = “2 and . Leotd , 
RS ds A 2k, ds 7 R, Rk, 
This varies between 0 when ¢ 37 and about 1/70 at the spherical 


junction. The solutions are therefore 
X = Ad-ieltisd J Aas, (16) 
Since expressions for s involve elliptic integrals, it will be necessary to form 
an interpolation formula for A as a function of s. 
Taking the equation of the hyperbola in the parametric form x = a cosec y, 
y = beot x, with k = a/b = tan#, we have 





x x 
* Teos?y l ay a ° | 1—sin28 sin? |! 
S a : X : cot®|~ dy ae | | ———. x! dy 
J [sinty © sin*ty sin? . sin*y ‘ 
Aor bar 
a eae ms , 
79 cot x(1—sin?? sin? y)* + cos?F | (1—sin? sin?y)-! dx 
sin 7 


har 


x 
| (1—sin?# sin?y)* dy |, tt 
dar 


on integrating by parts. 
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Thus, numerically, for the hyperbolic are between the cylindrical and 
spherical junctions, s varies between 0 and 20-5. 
Further, 
3(1—v*)|* (3(1—v?) |! 2 — 
—- | {(k?+-1)sin?d—1 
/ (Ry t) , (kat) 


Substituting for k, d in terms of # and x, we have 


tand = dy/dx 1 /(k cos x), 


[3(1—v2)} | 


(at) N (k? cot? -+-cosec*x) 54 (k? cot®x-4 cosec?y) | 


(18) 
Timoshenko (1, p. 482) suggests an interpolation formula for A in terms 


of s of the form 
Mt 


A (19) 

q Ss 
Since at the cylindrical junction s = 0, ¢ = 42, A 1/5, we must have 
mig = 1/5; and at the spherical junction s = 20-5, R, = the spherical 
radius 156 in. Hence g = 35 in. and m = 7. It can be verified numerically 
that in the interval 0 8 20-5 in. the interpolation formula (19) then 


wrees with the expression (18) to within 1 to 2 per cent. 
The solutions (15) for U or V then take the form 


A const x (¢+<s)*(q+s8)'tto™, 
or, With constants of integration A, a general solution is 
X (q+s 1,(q+-s)” 1.(q¢+8) mi | 
(q+s)-™*4| A,(q+s)'"+A,(q+38) aa (20) 
Writing 8 = mlog(q+-s)/q, dB/ds = m/(q+s) = A, an alternative form of 


the solution is 


X t -| e?( B, sin B+ B, cos B)+¢ 8( B, sin B B, cos 8). (21) 

\ q 
It is convenient here to take X U., and to derive V from (2) in the form 
R, 2U sais 


Et ds? ~ 
For the cylindrical and spherical regions, ?, is constant, and thus solutions 
take the form (15a). Considering the hyperboloid only as joining sphere 


nd cylinder, decay functions only need be taken for these two regions. 


3. Continuity components of stresses 
The solutions of (2), (3) obtained for the three regions considered, 
cylindrical, hyperbolic, and spherical, involve 2+4+-2 8 constants of 


ntegration, which have to be determined from boundary conditions at the 
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two junctions of the hyperboloid. For the shell of Fig. 1, which has a 
constant thickness, the following four quantities must be continuous at 
each boundary: 

(1 
(2 
(: 
( 


2) meridional rotation, V, 
) 
4) radial shearing force, Q. 


) hoop force, ogt, 
) 


3) hoop bending moment, ,, 


There are thus eight conditions for eight constants. The conditions 
of continuity of o,, and of the meridional bending moment being then 
automatically fulfilled. For vessels with t/R, < about 0-06, the meridional 
stresses do not approach critical values and can be neglected, so that we 
can now calculate the total hoop stress at the inner and outer surfaces of the 
hyperboloid. 

The above continuity conditions apply to the total values of the forces, 
which generally consist of two components. We will give the name ‘con- 
tinuity stresses’ to the component given directly by applying Hooke’s law 
for all parts of the shell where R, and R, have continuous derivatives. 
Where R, is discontinuous at a boundary, the material unity of the shell 
prevents a discontinuity in the expansion; this is achieved by the action of 
boundary forces acting in addition to, but depending on, the internal 
pressure. These cause additional ‘discontinuity’ stresses. 

To derive the continuity components, Hooke’s law gives the strains 
(1, p. 371) as 


€, = (o,—vo,)/E o,|1 v(2— R, R,)| E, meridional, 
and €, = (o,—vo,)/E = o,{2—v—R,/R,]|/E, circumferential. 


Then V is given in terms of the displacements u, w by the equation 








R,V = u+dw/dd, (23) 
where du/db—ucotd = R,«,— Ry, €, 
and w = ucot d— Ry €. 
Hence 
V = Rie cot $— © (Rae) 
O71 i. \ R, -_ d 9 
= —(—?__ ]|(—? 1+ 2v]cot 6 —— (Ry «,). (24) 
ag \z v}eo b— G7, Mee) 


For the hyperboloid, applying (1) and (9), the continuity component of the 
meridional rotation reduces to 


V, = ae of. (a s)eot d. 25) 
E\v? Un 
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The circumferential bending moment is 


M, p(ere V+y a (26) 


For the hyperboloid, 


hai \/1 , 
V cotd Fits -— 3}cot*d, 
R, 2Et\ * v2) \v 


and, using the expressions (5) for d R,/ds, and for v, we find that 


dV. (kK?+-1)(_ | 4/1 1/1 \) . 
: ) 08" 3 3) — —s }. —| - 3}}, 2 
ds * 2Et \' | (oa \(;2 ) “ls | 7 mee ) at 


Dp { (k?+1)/1 — .)\ — i. v{3 
| Vv o)+{k* Vv —vod})——|—— ‘ 
My 2 Eit\ v2 & ( a | ait Ie } sis y] 


(27a) 





Finally, the continuity component of the radial shearing force Q is given 
by (3), in which the left-hand side is taken as R,d?V/ds?. This reduces to 


» (k* sin 24] 3 /¢ (24-2 
R,Q _ Dp (kP+1 , in 4 ( , i) _ ih +2 pl. 
2 2 p2\ p2 


28 
2Et v \ 





The continuity components needed for the boundary conditions are now 
given by (1), (24), (27), and (28). 

The corresponding discontinuity components are determined by the 
solution (21) of the differential equation for Up, = R,Qp. The hoop force 
1, pp. 455-6) has discontinuity component —dU’/ds, corresponding to a 
component —t~!dU),/ds for the hoop stress. 

The discontinuity component of V is given by (22), and that of the hoop 
bending moment /, is found from V by (26): 


f 2/7 
o , [ R, ¢ ») Leotd 


: dU, 
Et ds\ ~ ds 


ds? 


Msn (29) 








4, Boundary conditions and numerical solution 
To apply the boundary conditions to the sum of the pairs of components, 
we require expressions for the corresponding quantities in the cylindrical 
ind spherical parts of the shell. In each of these parts A is constant: call the 
values A, for the cylinder, A, for the sphere. Using only solutions eventually 
decaying away from the junctions, expressions for U are 
1 


U, o—Adls—s (Ag etas 8’) | A,e iAd(s 8)) | 


Ty ‘ +s | ¢ id\s_| Aye iAys) 
(30) 


where 


A, m/q \[ 3(1 v*) |(at) + 
Ag = 4/[3(1—v?)]/,/(Ryt) = m/(q+s’) = m/qo, say, 


total length of the hyperbolic are. 
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For differentiation it is convenient to shift the origin of s and write s 
instead of (¢+s). Then s = q, q+s’ = q, at the junctions with cylinder 
and sphere respectively. From (20), with change of notation of the constants 


ai 


D C, gmA+D+44 0, g-m+d+4 4 c em(1—i)+4_1 0 8 m(1—i)+4 


and we shall write the terms in the first line as U,,, those in the second 
as U,,. Then 


d? _. (2am?2—1}) _, ee x: (20m?+-4) _, ; 
— Ups = + —— Ups: = Un = — ——- Up, (31) 
ds? s" ds? 3 





and in the expressions on the right-hand sides of these equations we can 
neglect the } in comparison with m?. Again, R, = as?/q?, so that 

277. D792 
2 Up ima 


R, Pit aeteed | 


2 2 D2 Up): (31a) 
ds q° 


Thus for the boundary conditions at the junction of cylinder and hyper- 
bola, the continuity of U and V gives 


C.+-C_.—C,-—C", A,+A, | 
2am?a , -~w a -" - an 
—,— (C,4+C_.—C,—C_,) = adj 2i(A,—Ay) | 

q~ 


where Cy, = C,qm+0+t, O_, = C_,q-™4+i+4, and so on. Since A, = miq, 


(32) 


the second of these equations is 


4+ < 


2 2 1 


1 = A,—A,g. (33) 


In these two equations we have used the facts that the continuity com- 
ponents (25) and (28) of V and R, Q vanish when ¢ = 47. 

Using (1), at this junction o,, = pa(2+-k?)/2t for the hyperboloid, and 
pa/t for the cylinder. Then for the hoop force, dU/ds, the continuity 
condition is 

(m(1+7)+-$}0,+{—m(1+i)+F0"_,4- {44 m(1—ic}4 

(4 —m(1 iC: 1 Lpak?q+-m{A,(1+7)+A,(1—7)}. (34) 

In the expression (26) for J/,, which applies to both continuity and 
discontinuity components, the term DV cot d/R, has been made con- 
tinuous by imposing the condition (33). Thus it only remains to satisfy 
the continuity of —Dv dV /ds, i.e. of the derivative of R,d?U /ds?, namely, 
using (31a), (27), 


. 2 . " 
($+m(1+7)'C5,+{4—m(1+i)}C_, 


| m(1 iC; i} m(1 iC" 1 


, 3.2 
= m{A,(1+i)—A,(1—iy} + "PE (a2 + 1y(k2—3). (8) 


4am? 








STRE‘ 


By 
(33) it 
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By adding (34) and (35) and then eliminating A, by the sum of (32) and 


(23) it follows that 


. 3 ].2 
C,4+{1—4m(1+i)}C"_, = kpak’g+ “P24 a k2-+-1)(k2?—3) 
q r " 4am- (36) 
° 00 
. 3].2 
C; ‘| tom(1 1c" ' } pak*q = | 4 k24 1)(k? —3) 
™ 4am- 


The conditions at the spherical junction can be similarly treated: it is 


convenient to write Cy = C,q7+o+4t, C*, 


= C_.qq ™+%+*, etc., and to 
denote by (R,Q.)o, Vy, d. the values of the unsuffixed quantities at the 
junction. 

The four conditions are therefore: 


for 1 CL4+-C.g+O1+07.1+( Be Qe = Aste: (37) 





1g | l 
MM | R, cot o,( - 
i \v. 


21m*a 


+1)(- 3) A,—A,; (38) 
* 


$+ m(1—i)}C}+{3—m(1—i)}Ct, 


non 
porn 


m(1+2)$C5+-44—m(1+2)'C*_,- 


i 
~ « 4 


spR,{2 \d0 kpqy R.—m{A,(1+7)+A,(1—2)}; (39) 


: d d?U\ : 
wd 107 | the f m | R, . of M,! x qo 2im2a: 
| ds ds? “| 
m-+-im)C,+(3—m—im)¢ 3 ($4+-m im)Cy (3 m-+-im)C” 1 
Etq?q,,| dv. ; . 
t fo) m{A,(1+i)—A,(1—i)}. (40) 
2im?a | ds |, 








Adding and subtracting as before, and using R, = aq?/q?, we get 
dV. 
ds 


1+4m(1+i):0%7+C", Ld, R,( 1 I 490 RY 


v5) 2m*a 


9 








(i —1) 4 _[ EtV],—(1+i)m(R,Q,)o. (41) 


2ma 


lhe conjugate complex equation holds for C; and (” 1: 


The numerical values of the constants can thus be found easily from 
36) and (41), in the form of numerical factors multiplied by p, the internal 
pressure. 

The calculations for the example chosen of a hyperboloid with a ratio 
a/b = V10 led to stresses exceeding 8 tons sq. in. on the outer surface of 
the hyperboloid. In order to limit the maximum stress to 8 tons/sq. in., 
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calculations were carried out for a hyperboloid with a lower a/b ratio. In 
order to use Legendre’s tables of elliptic integrals, it was convenient to 
choose # = 71° in (17), so that k = tand = 2-9042 = a/b. This hyper- 
boloid has a minimum meridional radius R, = —7 in. next to the cylindri- 
cal junction, and the relation (19) has m = 8, q = 40. The whole hyperbolic 
are is of length s’ = 24-5 in. so that g, = 64-5, and on evaluating the 
constants C,, C_,, Cy, C_, the expression for Uj, may be written 


m—4 


s 


8 
qo 


where 8B = mlog(s/q) and the constants B have the values: 


m+4 
(B,sinB+B,cosB), (42) 


B, = —270-74p, B, = +373-8p, 
B, 537-4p, B, = —I1I]p. 


Our particular interest is in the circumferential direct and bending 
stresses. Along the hyperbolic arc, which is the critical part of the shell, 
the direct hoop stress has the two components 


P(e = “3 _ pka {1+(k?+-1)sin*¢} 
ot R, 


2t {—1+ (k2+ l)sin®d}’ 





(43) 


and —(1/t)dU/ds. To obtain this, the length s of the are has to be deter- 
mined by (17) for each value of ¢. Then the second component of the direct 
hoop stress is 


1 « 
t ds : t j 


m—4 





Side as a ae ec —— 
(373-8 cos 8 — 270-7 sin B) (373-8 sin B+ 270-7 eos | 4 
8 Ss 


1/s\™+4#[m+4,_,- ; ne | 
+. “ (537-4 cos B— 11 sin 8) +- —(537-4sinB+1lcosf)|. (44) 
t\Go 8 ” 

The bending hoop stress is + 6/¢? times the bending moment M,,, whose 
two components are given by (26) and (29). Using (27a) and (31a) the 
total bending moment is 


Dp (v(k?+-1)/1 \ , [po =\E . « v (3 | Dp 2m? 
_ a a eS) re a ae hee to =e X 


m—t 
x {s50-6(4) [{q? cot 6—(m— })vas\cos(B —54° 5’) —mvas sin(B—54° 5’)]- 
Ss ; 


g\m+s_ ; . , 
537-5(= [{q2 cot d+ (m-+ })vas}cos(B+ 1° 11’)—mvas sin(B+ 1° 11’')] 
¢ 
) 


1 


| 
(45) 
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Fig. 3 shows the stress values for the outer and inner surfaces. Stresses 
determined by the support of the vessel have not been included, but the 
bending moments produced in the hyperboloid by the weight of the upper 
cylinder are opposite to those due to the internal pressure. Consequently, 
the stresses on the outer surface would be lower than those shown. 





Fic. 3. 


These stresses may be compared with the direct hoop stress on the 


membrane theory for the 10-ft. cylinder, and for a circular torus at the 


junction, of radii R, = 60 in., R, 6 in. (Fig. 4). The values are, for 
, ; ' R, R, 
the cylinder o, = pR,/t, and for the torus o, = p- +2 2), with a mean 
' a . t R, 
value 3-5pR,/t = 6 tons for p 45 lb./sq. in. 
~< O aM | 
| 
\| 
$ 
‘ Spherical } 
portion i, 
, ‘ J 
#) "il 


Fic. 4. Original design of bunged head. 


The direct hoop stress increases up to the junction with the sphere, where, 
lor the torus, R, = 156, R, 6, co, = 64 tons/sq. in., and for the sphere, 
To 2} tons/sq. in. 

[he mean of these two, 33 tons/sq. in., considerably exceeds the total 


hoop stress at any point for the hyperboloid. 
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With regard to the accuracy of the method, the continuity components 
are exact consequences of Hooke’s law. The discontinuity components, 
derived from the differential equation discussed, involve the approxima- 
tions on which the accuracy of the calculated stresses depends. This could 
be estimated from the approximation to the operator L in the equation (13) 
and the interpolation formula (19). In the critical part of the hyper- 
boloid, the discontinuity component of the total hoop stress does not exceed 
10 per cent. of the continuity component, so that the calculated total hoop 
stress remains correct to within a small percentage. At both ends of the 
hyperbolic arc the discontinuity component forms a larger fraction of the 
total hoop stress; but this does not matter much, as the total stress near 
the junction falls essentially below the critical limit. The maximum 
accuracy is therefore achieved in the middle part of the are where the 
maximum stress occurs. 
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SUMMARY 


The general theory of the deformation of a bimetallic disk with a small initial 


irvature due to temperature changes, is developed and applied to the case of a 
lisk initially in the form of a spherical cap. It is shown that if the central rise of the 


lisk is initially greater than about twice the total thickness, a position of instability 
will occur on both heating and cooling. An expression is derived for the mean 
mperature which a disk of this type will maintain when used as the control 


element of a thermostat, and results are obtained for the temperature variation 


ibout this mean value. 


Part I 

1, Introduction 
TuE control element of a thermostat sometimes consists of a bimetallic 
lisk which has a small curvature. If this initial curvature is large enough, 
the disk, on being heated, reaches an unstable equilibrium position and 
icks through to a stable one. On subsequent cooling a second position of 
nstability is reached and the disk clicks back to the original side. This 
feature means that a very positive switching action can be obtained by the 
ise of a bimetallic disk. However, the temperatures at which the two in- 
stabilities occur are different from each other so that this snap action can 
be obtained only at the expense of temperature hysteresis. 

The problem of how a bimetallic disk of this type deforms when it 
indergoes a change of temperature is considered in this paper, and a 
general theory developed. This is applied to the particular case of a disk 
pressed into a spherical form, and it is shown that, in order to ensure a snap 
ction, the initial central rise of the disk must be greater than approxi- 
mately twice the total thickness. This result is quite independent of the 
radius of the disk and very nearly independent of the particular combination 
imetals enrployed. Secondly, a formula is obtained for the mean tempera- 


ture which the thermostat will maintain. This is proportional to the 
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initial central rise and the thickness, and inversely proportional to the 
square of the radius and the difference between the coefficients of expansion, 
If the two metals have approximately equal thicknesses, the mean tem- 
perature is almost independent of their stiffnesses. 

The differential equations involved in this problem are essentially 
non-linear, and in order to obtain the temperatures at which instability 
occurs, solutions were obtained with the aid of a differential analyser 
designed by the Mathematical Instruments Section of the Common- 
wealth Scientific and Industrial Research Organization. These results 
enable the amount of temperature hysteresis to be determined in any 
given case. 

In the general theory, the only assumptions made, apart from those 
involved in large deflexion theory of flat plates, are that a temperature 
exists at which the disk is free from internal stress and that the disk is 
axially symmetrical. In the particular application of the theory it is also 
assumed that the Poisson’s ratios of the two metals are the same. 

After the work described in this paper had been completed it was found 
that Panov (1) had considered the same problem in 1947 and had derived 
differential equations of a type similar to those obtained here. The differ- 
ences are due to the fact that Panov made some additional simplifying 
assumptions which do not appear to be necessary. He showed, by an 
approximate treatment of his equations, that it is possible for unstable 
equilibrium positions to develop, but did not pursue the matter further. 


2. Basic theory 
Suppose that the disk, of radius a, is composed of two materials whose 
thicknesses are t, and ¢,, whose Young’s moduli are HZ, and F£,, and Poisson’s 
ratios o, and o,, and having coefficients of linear expansion a, and ag, as 
shown in Fig. 1. The total thickness (¢,-++/,) will be denoted by ¢. It will be 
assumed that at a certain temperature, which for convenience will be 
defined as the zero of the temperature scale, the disk is free from internal 
stress and slightly concave. Let w,(r) then represent the distance between 
a point in the weld surface at radius r and the tangent plane to the weld 
surface at the centre of the disk. Thus w,(r) specifies the initial shape ol 
the disk and is to be regarded as known. Throughout the analysis we shall 
assume that ¢ <a and that the range for w, is of the same order as f. 

If the temperature is increased by an amount 0, the curvature of the 
disk will either decrease or increase, depending on whether «, > or < 2%. 
During this change of shape, let w be the radial movement of a point in 
the weld surface at radius r, and let w be the final distance of the point 
from the tangent plane to the weld surface at the centre of the disk. 
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Then, on the usual assumptions for large deflexion theory of thin plates 
2). the radial and circumferential strains at a radius r and distance z from 
the weld surface (Fig. 1) are given by 


du ea : ] dw, - du’ 
3 dr ' (ar) “2\ dr)” dr | 


: l 
u Zz du’ | ( 
- r r dr / 
where Ww W— Ws. (2) 
MATERIAL 2 AXIS 








/ 
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j/ 
/ 
MATERIAL I. 


Fic. 1. 


Equations (1) give the total strains. Due to the temperature change 
alone, both the radial and circumferential strains are equal to a#. The 
strains due to the stress change are therefore given by 
e = ¢ xO. €g = eg—a, (3) 


? 


ind the radial and hoop stresses are related to the strains by the equations 
‘ 
5 (€9 T oe,). (4) 


On substituting (1) and (3) these become 





EK [du . ‘{dw\? 1/dw,\? u Ex 
a } | | o — 
] o~ | dr 2 dr 2 dr , r ] Co 
Ez d2w’ ' ao dw’ 
l—o? | dr‘ rdr|’ 





f. EB _ | u ‘ du f oT (a2) | : Kad 


l1—o?| 7 dr  2\d1 2\ dr, l—o 
Ez {ldw’ d*w" . 
; +o——}- (5) 
l—o*|r dr dr 





Che terms H, o, and «in these expressions take suffixes 1 and 2 according as 
218 positive or negative. 


5092.2 
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Let 7. and 7) be the radial and circumferential tensions, and M, and M, 
the radial and circumferential bending moments per unit length, so that 


ty ty 

T.= | f, dz, Ts | fo dz, (6 
4 “, 
ty ty 

M, | 2f, dz, Mp | 2fg dz. (7 
‘, ‘, 


We now define the following quantities for the disk: 


B- i, ty R o, E,t, 1% Et, 


' 
l ot ] o3 l—o;_-—sési - o3 


: a 1 /o, E, t? o, Et 
: 2\ 1 ot l a} 


D’ ] (" E, 4 , Bs E, :) 


9 1 € -~ | 9 
o; l—c 3\ l—oj l—o3 


Q 
l—o, l—o, Zz 


P E,ot, , B, Ya ty sq ~uG Hi, =) 


| Oo; ] O» 


Then, on substituting equations (5) into (6) we obtain 


T, Ble fe sie ee, ic 5 ctl 


dr ' 2\dr dr r dr” r dr’ 


T B Ui Bo (i) . (2) | Pe_C Ldw or Pw wi 


Yr dr dr 2\ dr r dr dr2 * 


Also, substitution of (5) into (7) and the use of (9) and (10) yields the ex- 


pressions 


C PC) C?\ dw’ CO'’\ 1 dw’ 
M. 7 p—1| Pp Dp’ nig 1] 
a (2 B ( s) dr? ( B F dr? 


C PC C2\ 1 du’ CC'"* d2w’ 
M, = —TM—(Q- §@—|D —{D’ — » (2 
2s” ( B ) ( 5) r dr ( B dr 


Equations (9) and (10) can be solved for w and du/dr, and on eliminating 


u between them we obtain the equation 


d 
dr 


d { d?w’ Ldw’ dw\?  [dwy\* 9 
3C’— B’C ; — | 3(B?— B” —i—)} |. (& 
(J . 5 =) r dr | 29 (ce) ( z) | 


B\T.— (rT)| = Bn w. (7,)| 


dr 
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Now, in the absence of external forces, the equations of equilibrium of an 


element of the disk are 


d 
T, = —(rT') (14) 
ar 
d yy dw 
nd My = —(rM,)-+rT, + (15) 


ibstitution of (14) into (13) now gives 


Bir? 37 |Z; (Bt BC) 





d? d l a 


r i se 
dr2 dr r|dr 


du:,,\2 du\? 
1( B2— B’? ? ; 18 
ye—w[(T) (iT) | a8 


\lso. on substituting (11) and (12) into (15) and using (14) we obtain the 
- d? d | 
|D \|, — 
B dr* ° dr r 


Equations (16) and (17) are two simultaneous differential equations relating 


equation , 
du dw 

¥ , 17 
dr ' "dr ey) 





T. anddw dr. and it is now advantageous to transtorm the variables involved 


ynon-dimensional form. Thus let 








ldw fan C2)]4 
- ae “dh. 
r di a | b b” 
l du, | (2(bD—C*)]3 
> i. “ho; 
] d a- Bb, Bb’? 
r- {- 
eS Pes © 
Ww. Ss 
(0-3) (18) 
Equations (16) and (17) then reduce to 
a. 2 “—_—" 
7 5 (9) }> 7 sind dy) | db b> (19) 
} d* 
d t “[n'd—d)] = oy, (20) 


a ; 2( BC’ — B’'C)}* (21) 
(B2— B’)( BD—C?) 5s 

he function ¢s(7) now defines the distribution of radial tension in the disk, 

the final shape and do n) the given initial shape. 
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If the disk is not constrained by any external means, the radial tension 
T. and radial bending moment M, are both zero at the edge r = a. On using 
equations (2), (11), and (18) these boundary conditions become 


g=6 a y=! - 
2(B Ta B2— B’2 © 4 
a*(BQ— PC ) XBD- raa| 6 
BD'—CC’. .d 
~ Bp aay at . @ 
| BD—C? — dy (pb—dy) a n l. (23 


The two remaining boundary conditions necessary for obtaining definite 
solutions of (19) and (20) are simply that ¢ and ¢& should be finite at » = 0. 

It should be noticed that the temperature @ only occurs in boundary 
condition (23). Thus the problem has now been reduced to finding solutions 
of (19) and (20) which are finite at 7 0 and satisfy (22). The tempera- 
ture corresponding to each solution can then be obtained from (23). 

A considerable simplification can be effected if we assume that the 
Poisson’s ratio of both materials is the same, so that 

O71 Do 0, say, 


B' oB, ag aC, Dp’ oD. (24) 


Equation (21) then shows that y = 0, and the problem is further reduced 
to finding solutions of the equations 





1° s 33 d* : 
4 (nb) = 8-82, 4 n($—d0)] = op (25 
dy? dy? 
with %=0 at 7 - (26 
l 2d . 
and i sib +3 sas (d—d¢y) at 7 l, (27 
3a2( BQ— PC —o* i 
where A= Ba(BQ—PC) _—s me he (28 
l-+o 2(BD—C?)3 


The assumption o, =: ¢, = ais most unlikely to give rise to any appreciable 
errors in the solution and it will be adopted hereafter. 


3. Disk initially in the form of a spherical cap 

To proceed further it is necessary to specify some initial shape for the 
disk. We shall assume in all subsequent work that the disk has initially a 
spherical form, so that to our order of approximation 


r y 
We = ho- (29 
a 


where h, is the initial rise of the centre of the disk. 
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From the second of equations (18) we obtain, after using (24), 
1—o? 


2 Bh x! 30 
| 3(BD—C?) we) 





Thus in this case ¢, is constant and equations (25) to (28) become 


d2 
f - ( nus) de d*, 
dy? 
d* 
4 = { np) dibs, 
dy* 
us(1) U, 
] 2 
Po J $(1) ¢'(1)}, (31) 
dy l+o 


where 7 is a non-dimensional temperature defined by 


dO (32) 
ho 
ind d’ denotes dd dy. 
It is easy to see that the following are three particular solutions of (31): 


a) ad do, f 0. T 0, 
h) db 0. us 1b2(] n); T :. 
c) f, do, uh 0. T 2. 


[he solution (a) represents the initial stress-free condition. Solution (b) 
shows that at a temperature defined by 7 = 1 an equilibrium position is 
possible with the disk perfectly flat, the radial tension, which is propor 
tional to ys, being then negative (i.e. compressive) throughout the disk. 
[t will be seen later that this state of equilibrium may be an unstable one. 
rhe solution (c) shows that at a temperature defined by 7 = 2, the disk is 
n equilibrium when turned exactly ‘inside out’, the radial tension then 
being zero throughout. 

It should be realized that the solutions (b) and (c) are a consequence of 
the initial spherical form, since they do not, in general, represent solutions 


of the more ceneral equations 25 


It will also be seen that if 


d =, ! . 


4 
Le) 


1 solution of (31). then so is 
dh @ bu TY, T 1+-f. (33) 
this means that we need only consider solutions between the limits repre- 


sented by (a) and (b) above. The solutions between (b) and (c) then follow 
immediately. 
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Let us now consider how the disk behaves when almost flat, i.e. when 
¢ is everywhere small compared with ¢y. In this case we can neglect ¢? 
on the right-hand side of the first of equations (31). Integration, and 
substitution of the conditions that (1) 0 and that (0) is finite, then 


gives 2 ‘ 
: ub 142(1 n). (34 
The second of equations (31) becomes 
d? e 
4 (nd) 1de(1—n)d. (35 
dy? 
Equation (35) cannot be integrated to give a closed solution. Accordingly 


we take 9 : 
, : db e(1- A,97+4, 7° ae (36 


where ¢€ (< dp) is a constant equal to 4(0) and the coefficients a,, are to be 
determined. Substituting (36) into (35) and equating coefficients of 1 


on the two sides of the equation gives 


d2 2 
a, adi As 2 ay), 
64 192 
> , 45 = 
a a, »—a ) n i: oe (31 
v ‘ 2 n—] 
; 32n(n+1)° ' 


Thus, for a given value of ¢, the coefficients a, can be readily determined 


expression 


€ 
T ] (38 
of € 
7 


O. > 
4. > (1 Wace Ja, : 
fe T i 








po 


It is convenient to express ¢ in terms of the central rise h of the disk. Fron 





equations (18) we have : 
db dy 1 
h 0 | r 
- “i | db dy. (39 
ho : do. 
| dy dy 0 
0 
Using (36) this becomes 
h “a 
: LJ S ite (40 
A: Py komt n-- I 


Combining (38) and (40) we now obtain the expression 


' e 
i pal [] 

where k i > (t+ a, 
lto 





If we substitute (36) into the fourth of equations (31) we obtain the 
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It will be seen from (37) that as ¢y-—> 0, a, > O (n 1, 2,3,...). In this 
case therefore k > 1 and +> 1—(h/hy). Thus if the disk is initially very 


dat the variation of the central rise with temperature is linear, as is shown 


vy the broken line in Fig. 2 


14 


a 














EEE 
O , a2 = 
-2+t 








Fic. 3. 

Several values were taken for ¢, and in each case a sufficient number of 
were calculated from (37) to enable the coefficient / to be determined 

m (42). Two different values were assumed for o, namely 4 and }, 
nd the results are shown in Table 1 and plotted in Fig. 3. 
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TABLE | 


do ° 4 38 10 12 14 16 
| | | 
| | > | 
Rig = | I ‘800 203 ‘261 835 "503 | 2°175 
r=} 3: : 5 
ki I ‘796 170 *311 | 9go4 1°589 | 2°266 





It will be seen that there is a critical value of 4, at which the coefficient | 
in equation (41) becomes zero. This will be denoted by ¢,,. and its value is 
given by as ‘ae 
given by i.=808 ££ o=3 


8-76 if o (43) 


_—— 


In this case small variations of shape around the flat form can occur with- 
out any accompanying change of temperature. The flat form of equilibrium 
is then a neutral one, and we can expect the variation of h/h, with 7 to be 
as shown by the curve labelled ¢, = ¢,, in Fig. 2. The coefficient k is 
positive if dy < ¢,,, and negative if ¢, > ¢,,. Curves corresponding to 
these two possibilities are also sketched in Fig. 2. 

It is now evident that if d, < ¢,, all positions of equilibrium are stable 
ones. If dy > ¢,,, however, the portion of the curve between the points 
marked B and EF in Fig. 2 represents unstable equilibrium. Thus, on being 
heated,+ the disk would follow the stable curve A F B and then jump directly 
across to the point C. On further heating it would remain in stable equili- 
brium along the curve CD. On subsequent cooling, however, it would 
follow the curve DCE and then jump across to the point F,, after which it 
would again remain stable. 

If the disk is to be used as the control element of a thermostat, it is 
important to know three things, namely, the critical initial rise required 
to ensure a snap action, the mean temperature which the thermostat will 
maintain, and the sensitivity, which is governed by the amount of hysteresis 
or temperature variation. These three factors will now be considered 


separately. 


4. The critical initial rise 
We have already seen that if d, > ¢,, a position of instability will occur 
on both heating and cooling. Thus, in order to ensure a snap action we 


have, from equation (30), 


p. > Per | 2(BD—C*)]} 
a! 2 
0 - 7 : 
2B 1—o? 
+ For the purposes of this argument it is assumed that a, x». If the reverse is true the 


words ‘heating’ and ‘cooling’ should be interchanged. 
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On using equations (8) and (24) this gives 


2 db f1-4 246 4 | 294}! 
h, | Por l lm(4m + 6m + 4) l?m*| . (44) 
t | 24(1—o?) |? (1+-m)(1+-lm) 

where l= H,/E,, m = t,/t.. (45) 


Using the values for ¢,,. given in equations (43), the critical values of 
jt, corresponding to an equality sign in equation (44), are shown in 


Table 2 for a few values of / and m, and for values of o = } and }. 


TABLE 2 


Critical values of h,/t 


o rl 
4 I 2 
n 
: 4 71 1°85 1°87 
1 1°77 1°85 1°77 
I 2 1°67 1°55 I°7I 


From this it will be seen that a wide variation in stiffness and thickness 
ratios for the two materials gives a comparatively small variation in the 
ritical value of h/t. We may conclude that in order to ensure a snap action 
the initial rise must be greater than approximately twice the total thick- 
ness, irrespective of what combination of materials is used. 

It is worthy of note that the critical value of h,/t depends only on the 
stiffness and thickness ratios for the two metals, and not on their coefficients 


f{ expansion or the radius of the disk. 


5. The mean temperature of the thermostat 

Referring to Fig. 2, it will be seen that if the disk is used as a thermostat 
control it will maintain a range corresponding to 7, < 7 < 7,. We have 
shown previously, by equations (33), that the points B and E are sym- 
metrically disposed with respect to the point 7 1, which therefore corre- 


sponds to the mean temperature 6,,. Using equation (32) we have 


f 
Po 


A 
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and on substituting for A and ¢, from (28) and (30) and using equations (8) 
and (24) we obtain the expression 


6, = : hot , (46) 
Xy Ny a“ 
: 1 —Im?)? 
where kK 4(1 1c) ( nll ats (47) 
' 4]/m(1+-m)? 


and / and m are defined by equations (45). 


The values of K for o , and } and for a few values of / and m are shown 


in Table 3. 
TABLE 3 


Values of K 


In view of the fact that in most cases the two thicknesses ¢, and 1, will 
be very nearly equal (m = 1) it appears that the value of A is almost 
independent of the combination of materials being used. 


Part I] 


6. The sensitivity of the thermostat 


The upper and lower temperatures of operation of the thermostat corre 


») 


spond to 7, and 7,in Fig. 2. Thus, if the thermostat maintains the tempera- 


‘ 


ture within the limits (@,,+-60), we have 


o@ (7,—1)@ 


where 6,, is given by equation (46). 


m 
The value of 7,, depends on the value of d,, and in order to determine it 
we now require solutions of the two simultaneous non-linear differential 
equations (31). These were obtained with the aid of a differential analyser 
(3) designed and constructed by the Mathematical Instruments Section 
of the Commonwealth Scientific and Industrial Research Organization 
and located in the Department of Electrical Engineering at the University 
of Sydney. 
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First, the dependent variables ¢ and & were transformed to new ones, 


f and q, by the equations 


] nd, g ny. (48) 
INTEGRATORS 
I bi pit w x wm wi Bisiig 


~ oSfh33 + 


_ mer 
Shean 


7 






f 
gq VOutput ffOutput Seo 
Fic. 4. Schematic set-up of differential analyser. (Scale factors omitted.) 
Since ¢ and & are to be finite at 7 0, the differential equations (31) and 


the boundarv conditions become 


ds - 2 : 1 
4 J be J = 5 4 aj Ig (49) 
dn? n* dy" ”* 
f(O) g(0) 0. (50) 
g(1) Q, (51) 
df 
T l ; 2 J — (] a)f ‘ (52) 
Dy | a) dy a=3 


In describing the solution of these equations it is convenient to express 


them in the form 


ag 1 42 a | if df tg l 
t Po 7) ’ )- “4 ( 7)- 

dy j n- dy 4) 7° 
The differential analyser was set up as shown in schematic form in Fig. 4. 
The notation follows that standardized by Bush (4). Integrators I and I 
served to generate the reciprocal 1/7, which was used as the argument for 
Integrator III, the output of which was the integral 

Ff 
—- dy. 
‘a i 
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This integral was used in deriving not only the right-hand side of both 


equations but also the value of { 6d» as the output of Integrator VIII. 


0 
The value of ¢y was set with the aid of the gear-box G. Values of f and g 
were plotted against 7» in order to examine their values at 7 = 1. The 


value of df/dy at 7 
grator VII. 
The differential equations (49) contain singularities at 7 


1 was read from the displacement scale of Inte- 


0, and to 
overcome this difficulty a series solution was obtained in the neighbourhood 
of this point. The condition (50) is satisfied by taking 


f = nlbo+by n+b. 42+...) 
+-C, n?-+...), (53) 
where b, and c, are the values of ¢ and & at » = 0. 


On substituting (53) into the differential equations (49) and comparing 
coefficients of the various powers of » we obtain the relations 


: . 2 2 

8b, bolo SCy 5 — bs 
24b, = boc, +b, cy, 24c, — 2b), 

. , ’ - ¢ 2 
48b, by Co tb, C1 +b C9 48¢, — 2b, b,—b} 


80b, by C3 +5, Co +b¢, +55 Cy 80¢, — 2b)6,—2b, b,, ete. (54) 


It will be seen from these expressions that if b, and c, are given, all the 


remaining b, and c, are calculable and the functions f and g are therefore 
fixed. In general, however, the function g will not satisfy the condition 
(51) so that there must be a functional relation between 6, and cy. Con- 
sequently, for a given value of 4», a value was chosen for b, (or ¢,). A value 
for Cy (or b)) was then estimated and the first few b and ¢ coefficients deter- 
mined from (54). The values of f, g, df/dy and dg/dyn at » = 0-05 were then 
calculated from series (53), sufficient coefficients having been obtained 
to give four-figure accuracy. These were used as starting values for the 
differential analyser which solved the equations over the remainder of 
the range 0:05 < y < 1. Inevitably this first estimate would result in a 
non-zero value of g(1). Accordingly, further estimates were made and the 
procedure repeated until a combination of b, and cy was obtained which 
made g(1) = 0 to a sufficiently high order of accuracy. Values of f and 
1 
land the integral [ 4 dy were then obtained from the machine, 
0 
and the values of h/hy and 7 to which the particular solution corresponded 


df dy at 7 


were determined from equations (39) and (52). Two values of d, were 


considered, namely ¢, = 12 and 16, and in each case a series of solutions, 
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corresponding to different combinations of b, and cy, were obtained as 
above. The results are given in Table 4, and plotted in Fig. 5, for values 
ofo = 4 and }. The tangents shown in Fig. 5 are given by equation (41). 


TABLE 4 





1 

| | ban Se Tie: 
1) 1) fli) | 6 h/hg o=} | o=] 
Je 00 454 “92 077 1°065 1'O71 
71 723 | 59 | 133 106 | mls 
C 63 1°005 2°51 "209. ||_=—s«1*153 1-167 
° I°II7 1°176 | 3°59 *299 1*194 1213 
53 165 | 4°49 "374 I°210 | 1°232 
18 188 734 5°46 "455 1°189 1°214 
I 58 038 6°88 573 1°135 1°163 
04 7 1069 | = 812 677 1°055 1083 
4 | 05 7°72 5703 | 10°79 *399 “690 “715 
02 513 1°085 ‘66 ‘O41 1086 1°089 
24 37504 | 2°29 "143 1°287 1°301 
5062 | 3°84 "240 | 1°431 1°454 
5658 | 5°14 “321 1*502 1°532 
597 6-006 5°97 "373 1544 | 1°578 
P | C ° 5894 7°82 “489 1°575 1616 
I } 5°318 | 8-96 560 551 | 1°594 
#340 | 1016 | -635 1°50! 1°547 
2°624 11°50 *730 I°392 | 1°437 
+362 12°98 ‘Sir |) or*180 1°221 
7 5-06 | 14°64 ‘915 | 827 855 





The maximum temperatures 7,, as obtained from Fig. 5, are given in 
Table 5. 
TABLE 5 


Values of 7, 


12 10 
[°23 I*55 
in be 
Remembering that, by definition, when ¢ = ¢,,; 7, 1, the curves of 


7,— 1) against dy are shown in Fig. 6. It will be seen that 7, and therefore 
the temperature variation 50, increases with increasing ¢). For maximum 
sensitivity of the thermostat the value of (7,,—1) should be kept as small 
as possible, consistent with a snap action. This means that the value of 
$) should slightly exceed ¢,,,, or that h/t should be slightly greater than the 


right-hand side of the inequality (44). 
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Apart from giving the values recorded in Table 4, it was also arranged 
7 

for the differential analyser to draw curves of g and | dd» against 7. 
0 

From these, curves of y% and w/h, against radius were constructed for the 


two solutions marked with an asterisk in Table 4. These represent con- 




















\ 
\ INITIAL 
75} . SHAPE 
\ 
tol : 
(@) DISTRIBUTION OF RADIAL (b) SHAPE OF THE DISK 


TENSION 


Fic. 7. Conditions near instability. 


litions very near to instability, and are shown in Fig. 7. The distribution 
fradial tension is shown in Fig. 7 (a) whilst the shape of the disk is shown 
Fig. 7(b). For comparison the initial shape, which is the same for both 


ilues of dy when plotted as why, is also shown. 
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SUMMARY 

In this paper the approximations necessary for a mathematical treatment of the 
general problem are first discussed. The main approximation involved is in the 
relationship between foundation reaction and plate deflexion, and a brief survey is 
made of the relationships which have been suggested by previous authors. 

The case in which the loading is static and the foundation reaction is directly 
proportional to the local plate deflexion is discussed in detail for both infinite and 
semi-infinite rectangular plates. Fourier integrals are used to solve problems in which 
all plate-edges are simply supported. The difficulties which arise when other edge 
conditions exist are discussed, and a method indicated for dealing with clamped 
edges. The problem of a semi-infinite rectangular plate under a given distribution 
of shear and bending moment along its free boundary is solved. 

The problem of a uniformly travelling load on an infinite plate is discussed in detail. 
It is shown that there exists a certain critical velocity, beyond which stresses and 
deflexions become infinite. The ratios of the maximum deflexion and _ bending 
moment to their static values are expressed as functions of A, the ratio of the actual 
velocity to its critical value. It is found that while these deflexions and stresses are 
greater than their static values, the increase is small unless A approaches unity. 


1. Nomenclature 
w = deflexion of middle surface of plate. 
p = mass of plate per unit area. 
h = thickness of plate. 
ao = Poisson’s ratio. 
E = Young’s modulus. 
D = };,Eh3(1—o?)-1 = plate modulus. 
k = foundation reaction modulus. 
1 = characteristic length defined by /* = D/k. 
x,y = rectangular cartesian coordinates in plane of plate. 
€, 7» = reduced rectangular cartesian coordinates in plane of plate. 


M;, M, = plate bending moments per unit length. 
M;,, = plate torsional couple per unit length. 


Q:, Y, = plate shearing forces per unit length. 
q = applied surface load per unit area. 
p, = foundation pressure per unit area. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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MATHEMATICAL THEORY OF A LOADED ELASTIC PLATE 33 
2. Introduction and formulation of problem 

The problem of the deflexion of an elastically supported elastic plate was 
frst discussed by Winkler (1). Since then the introduction of concrete 
paving slabs for roads and aircraft runways has given the problem increasing 
practical importance. A large amount of published work has appeared, 
mainly concerned with providing the engineer with practical design 
formulae. 

To translate the physical problem of a concrete slab resting on an earth 

yndation into mathematical terms requires a certain amount of idealiza- 
tion. The uncertain and complicated elastic properties of soils, as well as 
the possibility of plastic deformation, sets a limit to the validity of any 


tl 
< mathematical solution of the problem, and in practice theoretical results 
ey is | must be compared with experimental ones. 

We shall assume throughout this paper that the slab may be considered 
sa thin, uniform, isotropic plate. When discussing dynamic problems we 
hich | shall neglect the rotatory inertia of the plate. With these assumptions the 
edgt leflexion w satisfies the general equation of motion 
ope 
— DV4w-++(p/g)- = q(x, y;t)—p,(a, y;t), (1) 

at 
ta 
al where g is the applied surface load, p, the foundation reaction, and V4 is 
\ding the operator ' ee 


etua | | a) - 
sa 6x?" éy? 

Mindlin (2) has recently extended the more exact plate theory of Reissner 
3) to obtain the equation of motion of a plate under dynamic loading when 
shear deflexion and rotatory inertia are considered, but no application of 
this theory to the present problem has yet been made. In view of the un- 
ertain properties of all real foundation materials, it seems unjustified to use 
the more exact equation. 

The problem is therefore to solve equation (1) under given boundary 
nditions ard surface loading, when the relation between p, and w has 


een given. 


3. The foundation pressure p, 

In all previous published work, even when dynamic problems have been 
treated, the relation between p, and w has been considered as a static one. 
By making this assumption we neglect the effect of foundation inertia and 
loundation stress waves. There is also the possibility that for some soils the 
pressure may depend on the rate of deformation. A rough correction may 
made for foundation inertia by using an apparent plave mass per unit 
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area greater than the actual one, in the same way that allowances are made 
for spring mass in vibration problems. 

The simplest relationship is that proposed by Winkler, in which the 
pressure p, is assumed proportional to the local plate deflexion w. Hogg (4) 
has treated the case of an infinite plate under certain types of loading when 
the foundation behaves as a semi-infinite elastic continuum, and also (5) 
when the foundation is an elastic solid of finite depth. It is doubtful, how- 
ever, whether these are more exact representations of a normal foundation 
than that of Winkler. Some general relationships between p, and w are 
discussed by Holl (6), who obtains formal solutions of the general problem 
when the plate is infinite in extent. 

[t is clear that Winkler assumes complete lateral discontinuity in the 
foundation, while Hogg assumes complete continuity. An ingenious repre- 
sentation has recently been made by Hetenyi (7) whereby an arbitrary 
degree of discontinuity may be allowed. While he applies this successfully 
to one-dimensional beam problems, it is doubtful whether the equivalent 
two-dimensional plate equations would be tractable except in special cases 

In the remainder of this paper we shall assume the Winkler relationship 
to apply. While this assumption is an empirical makeshift, it has the 
advantage of being the easiest to handle, and has in the past been found to 
give usable results. Writing p, = kw, where k is an experimentally deter- 
mined constant, equation (1) becomes 


> 


(Dv -h+-(p 9) =) ust) q(x, y;t), (2 
ar 
or, taking reduced coordinates € = x/l, 7 = y/l, where /t = D/k, we have 
(vs +-1+-(p/gk) yee n;t) (1/k)q(€, 3 ¢), (3 
\ C ey 
where V‘ is now the operator (- oe ; +) 
C C x4 Cc ~ 


4. Boundary conditions 

The time boundary conditions of equation (3) may be treated formalls 
by operational methods. Restricting our analysis to finite, infinite, or semi- 
infinite rectangular plates, and assuming all loads to be in the finite region 
of the plane, the spatial boundary conditions will generally belong to one ot 
more of the following classes. 
A. Simply supported edge 
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B. Built-in edge 


ow 
Here w 0. —= 0. 
on 
# Free edge 
O-w O*w o [o*w . C2w 
Here Lo —— 0, | - + (2—o)—~ 0. 
con* Cs“ on\ on* os" | 


D. Edge with pre scribed stress distribution 


Here a2. 2 
uM’ — M kr i. 2] 
on cs* 
Qn = Qu— Fut = WE (EL + (2-0) 
cs on\ on* Ccs*~ 


In all the above equations n, s are reduced coordinates, respectively 


normal and tangential to the straight boundary considered, corresponding 


< 


to either of the coordinates € or y. M’, and Q’ are given applied shear and 
/ n n > 


bending moment distributions. Q’, includes the effects of the edge torsional 


t 


couple V,., since it is not possible to prescribe the three quantities M,,, M,,, 


n 


(), independently. + 


Tiere is also the additional condition, for plates with boundaries at in- 
finity, that w and all its first four derivatives must vanish uniformly there. 


5. Solutions for various boundary conditions 
5.1. STATIC LOADING 
In this case « quation (3) reduces to 
(V4+ 1)w(€, n) (1/k)q(€, ), (4) 
subject to the boundary conditions given in section 4. We consider first the 
ises where all finite boundaries are simply supported. 
Boundary conditions of class A 


The problem of an infinite plate under arbitrary static loading has been 
solved independently by Wyman (8) and Hopkins (9). 


(4) into polar coordinates and uses Hankel transforms, the 


The latter trans- 
lorms equation 
ntinuity of w and its first three derivatives rendering differentiation under 
the integral sign legitimate. 
This problem may also be solved directly from equation (4) by the use 
fadouble Fourier transform. The method is similar to that used by Holl (6) 


lhe relati | iv be fou y standard work—-see, for instance, 8S. P. Timo- 


aw-Hill, 1940), pp. 85-91. 
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in treating dynamic problems, and the solution may be shown to be 





. . i(aé +8) > ; 
w(€,n) = (1/4n7k) | eee dad i | q(u, vje-Heu+B dudy, 
= te 


(5) 
where (u,v) is the applied load distribution in reduced coordinates. 

This solution may be used to generate the solution for the problem of a 
plate occupying the region 7 > 0, and simply supported on the boundary 
» = 0. By taking a loading q(u,v), such that q(u, —v) —q(u,v), the 
solution given by equation (5) satisfies the boundary conditions of class 
A on 7» = 0, and the solution is given by 


L x 


w(é, n) (1/7?k) [ [ 1 ee dadB im 


mn 0 


‘al u, vje-"sin Bu dudv. 





Applying this reflection principle again, we may obtain the solution for 
an infinite quadrant occupying the region € > 0, 7 > 0, and simply sup- 
ported on the boundaries € = 0, 7 = 0. The solution is 


x x 


* sin a€ sin By 
+ (a2 4 B)2 


The problem of a finite rectangular plate may be solved along similar lines 





dadB i j q(u, v)sin au sin Bv dudv. 


by the use of a double Fourier sine series, the boundary conditions rendering 
differentiation term by term legitimate at the boundary.+ 


Boundary conditions of other kinds 

If, however, the boundary conditions belong to the other groups, the 
integral and series representations of w will not in general converge uni- 
formly on the boundary. Straightforward differentiation under the integral 
sign, or term by term, is thus no longer possible. 

If the plate is finite, and the boundary conditions are confined to classes 
A and B, however, double Fourier series may be used. The method is an 
indirect one and is given by Hopkins (10). He shows that for the case of a 
plate with clamped or simply supported boundaries, the deflexion w may 
be expressed as a double series which is differentiable term by term at the 
boundary. This series is, however, not a Fourier series and the terms are not 
orthogonal. The coefficients cannot therefore be determined separately, but 
appear as the solution of an infinite set of linear equations. In practice, an 
approximate solution may be obtained by considering a finite number 0! 
terms. A similar method involving integral equations could be applied to 
infinite plates. 


+ A solution is given in 8S. Timoshenko’s Theory of Plates and Shells, p. 251. 
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In cases where the boundaries are of types C and D the above method 
becomes unsatisfactory. The reader is referred to the original paper for a 
discussion of this point. 

Another method of attack presents itself when one pair of opposite edges 
iseither simply supported or at infinity. In this case a single Fourier integral 
series may be used in the direction normal to these edges, and the partial 
differential equation reduced to an ordinary one. This may then be solved 
by normal methods, and the solution expressed as a single infinite series or 
integral. 

Semi-infinite plate under given boundary stresses 


Asan example of this method, we consider the case of a semi-infinite plate 


occupying the region 7 > 0, loaded by a given distribution of shear and 
bending moment along the boundary 7 = 0. 
Jquation (4) become ' 
Equ O S (V4 l)w 0. (6) 


subject to boundary conditions D. 


We assume that w may be expressed as a single Fourier integral 
w(é, 7) W(x, »)¢ iat diy. (7) 


Similarly we assume that the boundary couple 4, and the boundary shear 


Q) may be expressed by 


é kl] ( M'(aje-* dx, (Sa) 
- Q = kl | Q'(a)e-*€ da, (9a) 


respectively, where M'(«), Q’(a) are given by the reciprocal relationships 


x 


M'( xy) (] 2arkl*) M,, ei dé, (Sb) 

Q'(x) = (1/2kl) | Q, ei dé. (9b) 
Equation (6) becomes 

oft ew : 

20° —~+-(a#+1)W 0, 

c "4 c i 
nd the solution of this which tends to zero as n> +00 1s 

W(x, 7) Fy (a)et + Fy(a)e-, (10) 
where b = B+ iy, and B, y are positive numbers satisfying the equations 

y B= x", 2By i. 
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The boundary conditions D become 


Q'(a) -|F(Fa-@-o)atr)) , 
oN\ on” n=0 
whence M'(a) , 


Q' (x) = iL F,(«)b{b?+ (2—oc)a*}— F, 

which may be written, on substituting for b 
M'(«) = —(1—o)o*| F,(«)+-F,(«)]+i[K(«)— K(a)] 

Q'(a) = —{8+(1—o)o*y}[ (a) + F(a) ]—ify— (10) 0B} F(a) — F(a) 
(11) 


, b, 
(x) 


If we write 


G(a) = F,(«)+ F(a), 
H(x) = i[F,(«)—Fy(a)], 
equation (10) becomes 
W (a, n) = [G(a«)cos By + H(a)sin By ]e-””, 
and equations (11) may be solved for G(a), H(«), viz. 


G(x) — 2x0) —12y*— (10) ah M"(o) 








l land es 
H(x) = —2(1 o)aPyQ"( x) + {1- : ee o)a2y?} M (a). 
l L4(1—o)a2 y _ o)at 


A formal solution of the problem is Poe il given by 





fT G(ax)cos Bn+H(«)sin By le-”” da. 


If the edge is a free boundary, a similar method may be used. The case 
of a concentrated load at a free edge has been solved in this manner bj 
Westergaard (11), who also gives approximate formulae for practical 
application. 


5.2. DYNAMIC LOADING 

elatively little work has been done on problems of dynamic loading. 
A formal extension of the theory to cover suddenly applied loads is straight- 
forward, however, as the time variable in equation (2) may be removed by 
operational methods. Holl (6) has given general solutions for both finite and 
infinite plates under dynamic loading, when the boundary conditions belong 
to class A and the plate is initially unstressed and at rest. It should be noted, 
however, that his solution for the case of a finite plate is valid only in certain 
cases, since by taking axes in the centre of the plate and expanding w ina 
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double cosine series a symmetrical form of deflexion is assumed. It seems 
that in the general case it is better to take axes at one corner of the plate and 
expand win a double sine series. 

We shall treat here a particular case of dynamic loading which may be 


of some pract ical importance. 


Travelling load on an infinite plate 

We consider a load of magnitude P, uniformly distributed over a rectangle 
of sides 2c, 2d, moving with constant velocity v across the plate. 

If we take the x-axis in the direction of motion, and assume that the stress 
snd deflexion patterns set up move with the load, then by taking reduced 
coordinates € = (w—vt)/l, n = y/l, w becomes a function of €, 7 only, and 


equation (2) becomes 
(V4+1+ (pv®/l2kg) aye) (1/k)q(, 7), 
og™ 


or, putting A? = pv?/2l*kg, 





_- 
[vs + 1+2d?. z)(E, 7) (1/k)q(€, 7). (12) 
\ssuming that w and g may be represented by double Fourier transforms, 
we have . 
: w(E,n) = (1/27) [  [ W(a, Bet$+8” dad, (13) 
q(é.) = (1/2m) [ ( Q(a, B)e"F+P) dadB, (14a) 
where Q(x, B) (1/277) ( ( qé, ne i é+Bn) dédn. (14b) 


\ssuming further that all stresses tend to zero at infinity, equation 


12) transforms to 

(a?-+%)24+1—2)2a2} W(x, B) = (1/k)Q(a, B). (15) 
rhe applied load q(x, y;t) is given by 
¢< £—vt < C. 


q(x.yst) = P/4ed for | 


lransforming the coordinates this becomes 





(—e/l <  <ell, 
djl - 


Y < ) P ted tor 
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and equation (14 b) then gives 


P sin(ac/1)sin(Bd/1) 


na ieiani << ete me 





Hence, using equations (13), (15), (16), the formal solution is 


x 


sin(ac/1)sin(Bd /I)¢ a€+B0) dad 


vB{(a?-+ B2)?+ 1—2A202 





w(E, n) (P/472cdk) | 


. . 
x - x 





or, in real form, 
oO « 


* ? sin(ac )sin(Bd I)cos a€ cos By 





vB{(a2-+-?)?+ 1—2d20?} 





w(E, 1) = (P/x?cdk) | 


0 O 





dad. (17) 


As c and d tend to zero, i.e. for the case of a concentrated point load, we 
obtain 


w(é,n) = (P/n20k) | | eee ante 


0 O 


(18 


The deflexion w will be finite at all points provided the expression 
(a? B)? | J] — 9) 2,2 
has no real positive roots. This will be so provided that A < 1, i.e. 
v? < 2kl*q/p. 

We thus see that there exists a certain critical velocity given by the equation 
v2 = 2kl?g/p. For velocities above this, our solution (17) becomes infinite. 

The physical significance of this velocity may be made apparent by con- 
sidering the propagation of straight-crested waves over the surface of the 
plate. 

If we consider a wave of the form 

Ww A sin 2z(x—vt)/a, 
where v is the wave velocity and a the wave-length, equation (2) gives 
v2 = (g/p){402D/a?+a?k/4n7\. 
Hence v? possesses a minimum value 2(kD)!qg/p for a wave-length 
a 27(D/k)*, 

or, from the definition of 1, v2 = 2k/*g/p, for a wave-length a = 2z7/l. Hence 
the quantity A is the ratio of the load velocity to this minimum wave 
velocity. The analogy with surface water waves is noticeable, the founda- 
tion pressure playing the part of surface tension in creating a minimum 
velocity. 


As far as aircraft runways are concerned, the value of this limiting velocity 
is well in excess of all surface speeds likely to be attained in the near future. 
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The following representative values are taken from a paper by Hopkins and 
Tones (12), for a concrete slab 1 ft. thick, 


k 750 lb./cu. in. 


BE = 6x 108 Ib./sq. in., o 4, h 12 in., 
The equation lA Eh? /{12(1—o?)k} 
then gives / 33-75 inches, corresponding to a value of 1,235 m.p.h. for »,. 


We now consider the effect of the velocity on the magnitude of the 


leflexion and bending moment under the load. We shall consider only the 


bending moment V,, as the effect on J, will be similar. 


sin(ac/l)sin(Bd/L) 


Writing O(a, f) 7 , equation (17) gives 
yc 
w(0.0) (P/z*k) | O(a, B) dad — (19) 
J J (a?+?)?+1— 2070? 
M,(0,0) = (PP/x2) | [ i eo (20) 
1 (x?+ B?)?+- 1— 2A?02 


0 0 


Asc 


Both these integrals converge to finite values. 


and d tend to zero, 


0,0) remains finite while J/.(0,0) does not. Note that ®(a,8)— 1/l? as 


ind d > (), 


The deflexion under a concentrated load is given by 


: dadB 


moO. QO) 


(his integration is best performed by changing to pol: 


3)-plane. Straightforward integration then gives 


w = {P/2nkl2(1+2))} [ (1—p2 sin?) 


where 2 2A?/(1-++-A?). 
W/W {2(2+-2 cos*v)! n\l K(v)— F(v,} 
where wy P/8kl?, sinv U 


(P/w?Pk) | | age pa 


ir coordinates in the 


dé, 


If wy is the deflexion under static load, we obtain 


177), 


and K, F are tabulated elliptic integrals, as 


given by Jahnke and Emde (13). A graph of the variation of w/w, with A 


is given in Fig 


1, and it will be noticed that although w/w, —> 0 as A> 1, 


tis not until A is quite close to unity that the effect becomes important. 


The increase in w with velocity may seem at first surprising. It must be 


borne in mind, however, that we have considered our plate as an undamped 
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elastic system. In most physical cases the effect of damping would be quite 
noticeable, and this would certainly:tend to reduce the deflexion. 

















08 5 02 4 06 08 10 
Fic. 1 


Since the bending moment under a concentrated load is infinite, we cannot 
simplify the integral in equation (20) by letting c and d tend to zero, as we 
did when considering the deflexion. If, however, we denote by M7? the value 
of M.(0, 0) for static loading, we have 


ang D(x. 2/ v2) 2 1 1p 
M,(0,0)— M3 = (2PPA}/nt) [ [ ,-R'oPlate’ tek) Cate) 
: J J f(a + B?)?+ 1 (a?-+ B?)?+ 1 — 20202) 
0 0 
and this integral remains convergent. Since ®(a, 8) < 1/l? when the loaded 


area is finite, equation (21) can be simplified to give an upper bound for the 
increase in M,, this bound being approached as the loaded area decreases. 
Thus 

x? (a? +-oB?) dadB 
{(x?+ B?)?+- 1H (a?+ B?)?+ 1—2A?%0?} 


(22) 


M;—M? < (2P)2/n?) } | 


This double integral might be reduced to elliptic form, but it is simpler to 
obtain the solution as a power series in A*. Expansion of the integral in 
equation (22) gives, after some manipulation, 


x 


M:—M? < > (PA*"/2n?) [ sin”-!¢6 dd [ [(1—a)cos?”+26+ ¢ cos?”6| dé. 


n= 0 0 
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M,—M? < P ¥ a,*, (23) 
" n=1 
2 9 1 = ' — | 
where a “ 2)(gn—1)! (n— 2) (24) 
24r(n- 1)(4nr—4)!n! 


The first six coefficients are as follows: 

















a, 0-05195, a, = 0-01354, 
As 0-02652. a, = 0-01089, 
As 0-01790., Ae 0-00912. 
1A 
' 
2F ' 
M,, / 
Ke : 
J 
Q t i oll. 1 4 
: 0:2 0-4 0-6 0-8 10 
A 
Fic. 2 
From equation (24) it may be shown that, for large values of n, 
a, ~ 1/(4av2n), 


ind the series is therefore convergent for A < 1. 

We now consider as an example the slab treated previously in considering 
the critical velocity. For the case of a total load of 50,000 lb., distributed 
veranarea of 500sq.in., Hopkins and Jones (12) give a value of 310 |b./sq. 
n., for the maximum bending stress under static loading. This corresponds 
toa bending moment MV? equal to 7,440 ft. lb./ft. [These values are given 
lor a circular ly loaded area, but it is not likely that the error will be large. | 

Equation (23) may then be written 

M./M? < 14-6-73 ¥ a, A". 
n=l 


\ graph of the upper bound of .V,/M? is given in Fig. 2. 


6. Conclusion 
In the case of static loading, two problems of practical importance remain 
to be solved. The first is the establishment of a relationship between the 
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slab deflexion and the foundation pressure which agrees more closely with 
experimental data. The second is the solution of problems, such as that of 
a loaded corner with free boundaries, where the boundary conditions intro- 
duce mathematical difficulties. 

The treatment of travelling loads in this paper must be regarded as very 
approximate, owing to the assumptions discussed in section 3. The results 
obtained, however, do show that for most purposes the effect of load 
velocity may be neglected. It should be pointed out, however, that the 
effect of foundation inertia will be to lower the calculated critical speed. 
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\ WEAK SPHERICAL SOURCE IN A ROTATING FLUID 
By K. STEWARTSON 


(Department of Mathematics, The University, Bristol) 


{Received 21 February 1952] 


SUMMARY 
The effect of placing a weak spherical source in a rotating fluid is discussed, and 
s shown that ultimately the whole output of the source is funnelled inside a 
nder circumscribing the sphere and having its axis parallel to the axis of rotation 


fluid. 


1. THE motion of a sphere along the axis of a uniformly rotating fluid, and 
the difference in the character of the ultimate motion of the fluid as com- 
pared with that caused by a sphere moving in an unrotated fluid, have been 
liscussed in an earlier paper (1). In this paper we propose to use the same 
techniques to determine the ultimate motion of the fluid when a weak 
spherical source is placed in the fluid. 

We shall show that, just as in the previous paper, the cylinder @ circum- 
scribing the sphere, and having its generators parallel to the axis of rotation, 
lays an important part in the ultimate motion. Inside @ the perturbation 

the velocity of the fluid due to the source is wholly parallel to the axis 
of rotation and has a singularity on @. Outside @ the perturbation consists 
fa transverse velocity about the axis of @, singular on @, and behaving like 

?as r (the distance from the axis of @) tends to infinity. A sketch of the 
streamlines is given in the diagram on p. 46. The funnelling of the output of 
the source parallel to the axis of rotation is in direct contrast to the spheri- 
cally symmetrical output when the angular velocity of the fluid is zero and 
provides the chief interest of this work. In addition to the steady picture 
lescribed above, the velocity components are oscillatory on the sphere and 

its axis, while on @ they increase without limit. 

The problem has also been discussed by Grace (2), who used expansions 
ff the velocity and pressure in powers of t, and obtained the oscillatory 
motion on the sphere and on the axis of @. He did not, however, obtain the 
general pattern of the ultimate flow. 

2. Let the undisturbed motion of the fluid be a uniform rotation with 
ingular velocity } about a fixed axis and let us take an orthogonal set of 
Cartesian axes Oxyz with origin O at the centre of the source and Oz 
parallel to the axis of rotation. Further, let the axes rotate about Oz with 
the same angular velocity as the undisturbed fluid and let the components 
of fluid velocity relative to these axes be (u,v, w). Then if the radius of the 
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46 K. STEWARTSON 


source is a, and it begins to emit fluid with constant velocity W whent = 9 
the boundary conditions are that 

zu+yo+zw = Wa when 2?+y?+22= a? (t > 0) 
u=v=w=0 when 2?+y?+22> a? (t= 0) (] 

and u,v,w>O0 as 2+ y?+22 > o9 for all t 
We now assume, with Grace, that it is legiti- 
mate to linearize the equations of motion. 4 

necessary condition for this is that 


{| | rT | W = o(a), (2 
ATALA| 4 


which provides not only an upper bound for 
W for fixed a, but also a lower bound on a for 
fixed W, so that, in particular, we cannot 
Lt | te determine the flow due to a point source by 





proceeding to the limit as a — 0, as was done 


ae, | 





by Grace. In spite of the restriction that IW 








must be small compared with a, we still 


N 
\ 

a 
© 


cannot justify linearization indefinitely be- 
cause it implies that the vorticity oscillates 









































ry infinitely. Comparison with experiment does 
NL not suggest that this is serious, however, pro- 
vided only that the predicted ultimate flow is 
physically possible. We write 
je an p — - 
NL LILY “alee 
a where R is the distance from the axis of rota- 
Ba a tion, p is the pressure, and p is the density; 
further, let a bar over a function denote its 
Laplace transform with respect to ¢, so that, 
AIT Hk for example, 
LLY 
1 
P = [ eP dt. (4 
Sketch of the streamlines in the _ 


perturbed motion. o ; an P 
Then in a similar way to (1), or applying 
the Laplace transform to the equations used by Grace, P satisfies 
2p 2p | e2 92p : 
o2/ é I 1--s? 07F o: (5 
9 1 ‘ 9 9 as b 
Cn" Oy a G2 


s oP 1 oeP l1 oP s oP 


stl dx s*+l1 dy s?+l dx s*+1 cy 
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n 1 aP 
ind Ww ais 
S C2 


Interms of P, the boundary conditions are that 





[ oP oP l+s? oP or oP Wa —" 
S| x y — 2 x — y— im s (1 + 8?) 
C2 cy a C2 cy Cx 8 
when 22+ y?-+-2? a and 
D P P 
Ca ¢ c P . ‘ ~ 
—s >O as 22+ y?+22 > o. (7) 
Ca cy Cz 


if the coordinates are now transformed to C, uw, @ (1, equation 2.16), where 


a a *9\1 9\1 ; 
pc, } — (1+ ¢*)*(1—p*)?, ¥ r cos 6, y rsin@, 
g ( g2) 
(8) 
he boundary conditions simplify to 
eP l oP Wa . 
when ¢ 8, (9) 
ray s- 1 ot S A 
and P—+>0O as C > CO, 


n which case a solution of (5) may be found which is a function of f only, 


iaW (1-+-s8?) iL -§ 
P log (5 }. (10) 
28s _=—"s 
Inverting the transformation to determine P, u, v, and w as functions of 
j,z, and t we find that 
Wa f 1-+8? C+1 
P | log |; Je ds. 
>) ~) ¢ “4 
mil A ats < t 
Wa | (Sa y)(1 s2)est ds 
/ — 
2m J 8 C*)(2a?(*— s?r2— 8722 1? +a?) 
(11) 
= | (sy—2)(1-+52)e% da 
v = - 
2m J 8 C*)(2a?l? —s?r?2— 222 1?-+-a?) 
Wea f° (1-+-.62)e% ds 
| { s : (Cc > 0). 
Dari 2a2l? — 272 — 5222 2 +a?) 
the general characteristics of the flow may now be found in the same way 


sin (1). We are really only interested in the ultimate flow, however, and 
‘0 we shall content ourselves by observing that the only contributions to 


the various functions which do not tend to zero as t > © are those from the 
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poles of the integrands; and, further, in general the only poles occur at s = ¢ 
We shall now obtain the contributions giving us the ultimate flow. 





(a) Outside @, r > a and (1, equation 4.2) f = (r?—a?)!/a when gs = 
Hence as s > 0, 
, W ya? = Waa? 
6 > aa oe easy (12 
r*(r*—a*)? r*(r*—a*)? 


while @ is bounded, so that ultimately the perturbed part of the velocity 

is in the @-direction and equal to 

Wa? r 

ra 

Similarly P > aW sin-!a/r as too. The motion outside the cylinder 

yr =a is therefore ultimately purely two-dimensional in character, con- 

sisting of a swirling motion about the axis of @, with a singularity on % 
and rapidly tending to zero as r > 0. 

(b) Inside @ but not on its axis, we have 0 < r < aand (1, equation 4.13 


9 


sz(a2—r?)-* near s = 0. In this region @ and @ are both bounded near 


wo 


s = 0, while sii > Wa(a?—r?)-4 and sP > 4naW. (14 


Hence ultimately the velocity has an axial component only, equal to 
Wa(a?—r?)-!, (15 
while the pressure term P tends to }zaW. 
In addition to this general pattern there are also special regions where 
poles occur at other points in the s-plane. 


(c) When r = a, and z + 0, we have (1, equation 4.6) Z (sz/a)! neat 
s = 0, and as before, when t is large, 
+ £ A ut i -(at\4 
umv u( )*. va —NV x{ —}°, and wz | -1*, (6 
“ \q7raz TAZ TZ 


so that the azimuth and axial components of velocity increase without 
limit on @. 


(d) When r = a, and z = 0, we have € = s and so 
W W . 
u (yt+2x), v (y—tx), and w= 0. (17 
a a 


The axial velocity is therefore zero at all times, the radial velocity is equal 
to W in agreement with the boundary condition, and the azimuth velocity 
is equal to — Wt and decreases uniformly with ¢. 
(e) Onthe axis of the cylinder wherer = Oand Z = sz/a. Hereu = v =! 
by symmetry and 
c+ia 
WW * (1+ 8?)a?e* ds w( z*7—q? “) 


w= 1 ———_ cos 


27t J = s(a?-+-8?2?) 
c—io 
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ind oscillates finitely about a mean value of W. Also 


at/z 
_ [ F—a®. 
- Wa —— sin a da. (19) 
3 ue" 
0 
On the sphere r = a where (€ = s and hence 
Wa - (y+as)e* ds Bf 2t yw. & 
u ated ee w( cos — 4 sin : 
27. a*s*--2* a a 2 a 
-[y zt y... & 
a WW ( cos an —i. 
\ a z a 
aw a*—2* zt 
Ww (1 cos ), (20) 
Z a“ a 
 f sing Wa , 
nd P av da +- ——|tcost sin ¢}. 
x ig 


« 
0 


[he velocity and pressure distributions in these last three regions were also 
found by Grace, but we can now see that they are special cases and give no 
suide as to the ultimate behaviour of the fluid in general. 

Like the problem of the uniform motion of a sphere along the axis of 
rotation, so here the streamlines in the ultimate flow are either straight lines 
reireles, and so there is a distinct possibility that the theoretical picture 
sin good agreement with experiment. If this is the case and coloured liquid 
were used in the source, the experiment would be very striking because only 
the region inside the critical cylinder would be coloured, the rest of the 
fluid being untinged. 

Note added in proof.) The boundary conditions given in (1) are not com- 
pletely satisfactory, because there is a discontinuity in the velocity distribu- 
tion when the source begins to emit fluid. As Morgan (3) points out, it is 
better to use the irrotational solution for a boundary condition when ¢t = 0. 
In this problem we ought to take 

U,V, W gradd, whent = 0, 
vhere d W (x?+- y?+-2?)-3. 
‘light changes in the detail of the solution would then have to be made 
but the results are unaffected. 


REFERENCES 
1. K. Srewartson, Proc. Camb. Phil. Soc. 48 (1952), 168. 
2.8. F. Grace, Proc. Roy. Soc. A, 105 (1924), 532. 
3. G. W. Moraan, ibid. 206 (1951), 108. 














THE FLOW OF FLUID ALONG CORNERS AND EBGES 


By L. SOWERBY 
(Department of Applied Mathematics, University of Liverpool) 
and J. C. COOKE (University of Malaya)t+ 
[Received 18 October 1951] 
SUMMARY 


The steady laminar flow of fluid along corners and edges has so far proved to | 


an intractable problem by using boundary-layer methods ; only one attempt, due t 


Carrier (1), is known, and this must be regarded as unsatisfactory in that certain of 


r4 


the equations of motion have been disregarded. It was suggested in Goldstein (2 


that qualitative results could be obtained by an hypothesis originally due to Ray. 


leigh (3), and it would seem that quantitative results may be obtainable by changing 
one constant. In Part I of this paper the modified hypothesis is tested on the floy 
along the outside of a circular cylinder, some results being obtained for compariso! 
by standard boundary-layer methods. The results suggest that some estimate of skit 
friction may be possible. In Part Il of the paper the problem discussed is the flow 
of fluid bounded by two semi-infinite intersecting flat plates which are both set 

motion suddenly with uniform velocity W parallel to their line of intersection (‘thi 
wedge problem’). The modified form of Rayleigh’s hypothesis is applied to certair 
of the results obtained to estimate the drag on a finite rectangular pipe of suitab 


length mounted in a steady stream. 
Part | 
1. Introduction 
IT was originally suggested by Rayleigh (3) that information on the stead, 
flow of fluid past a semi-infinite plane with its edge at right angles to th 


stream could be obtained by considering the problem of an infinite plane 


moved suddenly in its own plane with velocity W. It is easy to solve this 


problem in terms of the time ¢. The suggestion was that ¢ be replaced bj 
z/W, where z is the distance measured downstream from the leading edg¢ 
to obtain information on the flow in the steady motion problem. It is tru 
that there is thus obtained qualitative information which is reasonabl) 
accurate, but the drag coefficient for a plate of length / is found to b 
4/(7R)} 2-25676/ R!, where R W//v, instead of the correct value whic 
is 1-32824/R}. 

A change of constant in the transformation from ¢ to z could be made t 


yield the correct value for the drag coefficient in this case, and it seems 


possible that in more general cases one should be able to make the sam 
change in the constant. It is true that, to a first approximation, this chang 


In the first instance both authors submitted papers independently. However, with t! 


exception of Part I, which is due to Professor Cooke, a substantial amount of subject-matte! 


was found to be common to both papers, and it was decided therefore to combine ther 
into a single paper. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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an be made for the flow along a semi-infinite circular cylinder, as will be 
seen later 

In view of this, it seems worth while to consider Rayleigh’s problem for 
the wedge; that is, we consider the motion set up in a viscous incompressible 
fuid, which is bounded by two semi-infinite flat plates intersecting at an 
ngle 8 (0 < 8 < 27), when both plates are set in motion suddenly from 
rest with uniform velocity W parallel to their line of intersection. 

Certain special cases of this problem have been discussed already, begin- 
ning with Rayleigh’s solution (3) for the infinite flat plate. Howarth (4) has 


nsidered in detail the problem of the semi-infinite plate (8 27), and 
Sowerby (5) has considered solutions in the case B = z/n, n being an 
nteger |. Jaeger (6) gave results for the wedge problem, treated as a 


problem in heat conduction 

In this paper it is suggested that information about the related steady 
low problem might be obtained by replacing t by y*z/W, where y is to be 
letermined from the flat plate problem with its known solution. 

For clarity we state here that ‘length’ is to be taken to mean a distance 
easured parallel to the direction of motion, and ‘width’ to mean a distance 


easured perpendicular to the direction of motion. 


2. The equation of motion, and boundary conditions 

Using cylindrical coordinates 7, 6, z with the z-direction parallel to the 
lirection of motion, It is possible to find a solution of the Navier-Stokes 
equations in which the only component of velocity w is parallel to the 
direction, w being independent of z; these equations and the equation of 


ontinuity are satisfied if 


ot ‘N\ér? "+ or * r? 06? 
the pressure being constant everywhere. 
For the infinite circular cylinder of radius a moving parallel to its 
generators the boundary conditions are 


w=—O0 when t=0 (r>a) | 
when f 0 | 


‘being independent of 6. 

For the infinite wedge of angle B. placed with its corner parallel to the 
lirection of motion, the boundary conditions are 

w 0 when ft 0 (0<6< 8B) 

w Vi if @ 0. 60 8 

w>0 ar>o (0<60< 8) 


| (2.3) 


when f 0 | 
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These are problems in heat conduction, and their solutions may be found 
in Carslaw and Jaeger (7). 


3. The infinite circular cylinder of radius a 


The solution of this problem is given in reference (7), p. 281. For small 
values of ¢ the solution is 








w as » . (r—a)(vt)? . », , (9a®—2ar—Tr*)\it., 
= = =< @380(C)-L ——_ 2 sg erfo(7) 4. — 2 erfe(C)-+-..., 
VW r ($) 4 4atri st 32airi ($) 
- r—a 
where C = —, 
2(vt)! 
and i” erfe(x) a"-lerfc(£) dé (n = 1, 2....), 
” 
— 2 / ‘ 
Merfe(x) = erfe(x) = 1—erf(x) 1—— | exp(—€*) dé. 
ar 
0 
The skin friction 7 is given by 
ow W mvt)! vt at ' 
T me *} = | sen ——. + —~ (vt)? 4 ms | 
Or}raq (mvt) | 2a 4a?" 8a j 
We now apply Rayleigh’s hypothesis to this, but write ¢ == y*z/W instead 
of t = z/W. This ‘converts’ the problem into that of steady flow along the 
outside of a semi-infinite cylinder. If we consider a length / of this cylinder 
the total drag is , 
| 2nar dz. 
0 


Putting R W1/v and (vl/W)! ea we obtain for the total drag the value 





nal ll ll i 2 
yRt |i , —" a 2 


and if this is equated to }Cp 2zalpW*, where C, is the drag coefficient based 
on the total area, we obtain 


4 a yer et 
ee oth aks ies Lae he Pest _..). (3.1 
D yak 4 12 732? | 
Thus the first approximation for a short length of cylinder is C,, = 4/y(7R), 


which is the same as for a flat plate. Now it is known that the first approx- 
mation for the drag coefficient of a circular cylinder from the boundary-layer 


equations, provided the length is quite short, is the same as for a flat plate 
namely a/R! where a = 1-32824. Hence we put y = 4/az? = 1-699. 
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[It would be a good test of this hypothesis if it were found to give com- 
parable numerical results in the case of higher approximations, that is, for 
: longer length of cylinder. We proceed therefore to consider the steady 
dow in the boundary layer past a circular cylinder. The method has been 
ysed by Atkinson and Goldstein ((2), p. 304) in connexion with the flow 


nside a pipe. 


4, Boundary layer flow along the outside of a circular cylinder 


The equations for steady flow in cylindrical coordinates are 


ne : a 
ow ow l Op Ow low . dw 
u Law +v{—-+-—4 =I, 
or Oz p oz or ror oz 
° P ai 
ou ou lop (eu low, dw u 
u ~ 2 -vV naasye: Sipet - j— ——_— ——- —— J. 
' ‘ 9 ' ‘ ‘ 9 9 
cr C2 por or* ror Cz" i? 
ee Ow 
(ru) +— = UV. 
Yr o7 OZ 


Making the usual boundary layer assumptions, but not assuming the 
boundary-layer thickness to be necessarily small compared with r, we 
\btain the equations (for constant velocity outside the boundary layer) 


Cw Ow v ¢ Ow al C 
u w (r ; — (ru) -+—(rw) = 0. 
oO} Oz r oO? cor or C2 
Hence there exists a function % such that 
l lé 
Ww ad “= — ad 
ro? r Oz 
_ l ] { »2 
» 0,% . ra | 5 ? U 
Writing E | - “| Y1 =—1}, y= sy 
a*V 2\a? 2 
ind ob a? W {Ef,(n)—&fo(n) +8 f5(7)—...}, 


where W is the constant velocity outside the boundary layer and a is the 


radius of the cylinder, we have 
w= sWifi—&fet+efs—...}, (4.1) 


the dashes denoting differentiation with respect to ». The functions f 


satisfy the equations | " ' fot 0. (4.2) 
Sotho h—fefit2hfi = 4nfi+4fi (4.3) 


Istish —2f 3h 3f3 f 3 4nfo T 4f 2h fs tf, etc., (4.4) 
with f,(0) = f,(0) = 0; f\(00) = 2; f,(00) = 0 for s >1. Now the skin 
iriction is given by 


Cw LL W " wa orn 
r= nS) = BE p50)—2F40)+ €FH0)—..9, 


40& * 
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and the total skin friction for a length / is equal to 


Performing the integration and simplifying we obtain for the drag coefficient 


| a ” of” 
Up pis 10) bef (0) T 3€°f 3(0) sess 


X { ] - Bo yates ) 
|] — o(O) + “f a(O)—...}, 5 
R})\ 3,64 2 Sats al ) j (4 


where « and R have the same meanings as in the result (3.1), and 
x = f;(0) 1-32824, 


from the well-known solution of equation (4.2). 
The equation (4.3) has been solved numerically to a rough approximation 
. and it yields f3(0) <= —2-77. 

We now compare results (3.1) and (4.5). To make the leading coefficients 
agree we put y = 4/a7! as already explained. With this substitution we 
r test the next coefficient. If they were to agree we should have 


: my /4 £3(0)/2a, 
; which would give f5(0) 2. The agreement is not good, but at least it 
. is of the right sign and the right order of magnitude. Presumably in order 


to ‘force’ further agreement we should multiply the ¢ in equation (3.1) by 
a suitable scale factor, and test by considering further terms. This has not 
been done here, because there is not much point in trying to force agreement 
further. 


According to Seban and Bond,*+ f3(0) = —2-816 (instead of —2-77 as 
given here) and f3(0) 1-520. These would make the factor in brackets 


in equation (4.5) equal to 
1 + 1-060€—0-381e?4 
whilst that in equation (3.1) is 
1+ 0-753e—0-241e?+.... 
We proceed now to discuss the wedge problem, and finally to apply the 


hypothesis of Rayleigh to certain of the results. However, instead of putting 


t = z/W we shall write t = y*z/W where y = 1-699, thereby expecting to 
obtain results of fair accuracy. 


+ Journ. Aeron. Sci. 18 (1951), 671. 








5, Ra 

The 
the eq 
tively 
knowl 
sppea 
For tl 
may | 


Wr 


then 


where 


The } 
lefine 


ind t 
A, of 


carry 


is det 


Irom 





ficient 


atior 


cient 


on wi 


AST It 


orde T 


iS nol 


ment 














THE FLOW OF FLUID ALONG CORNERS AND EDGES 











Part I] 


5, Rayleigh’s problem for the wedge 


[he statement of this problem has been given in the introduction, and 
e equation of motion and boundary conditions in (2.1) and (2.3) respec 
ely. The problem is one of heat conduction; accordingly we can use 
nown solutions, and the most useful form of the solution we require here 
ppears to be that involving a Green’s function developed by Carslaw (8). 
For the theory of Green’s functions in heat-conduction problems reference 
LY be made to Carslaw and Jaeger (7) (chapter 13). 


Writing, for convenience, 


bh W uw. (5.1) 
B x 
n w*(r. 0 / W dé’ | yr’ dr’. (5.2) 
0 0 
Pre 
| { (r2+-7°*)) 
— xX] 
Savpt | byt | 
> r’ COS a’ | ] | - 
exp! Fi "i = da’. (5.3) 
2vt ; || etpa etVY—-O—) — pipa ei p+’) 


parameter p, which will be used frequently throughout this paper, is 


fined by the relation a 


p 7/B. (5.4) 
nd the contour A’ in the a plane consists of the two curved parts A, and 
{, of Fig. 1, which is taken from Carslaw’s paper. 

Carslaw obtains an expression for v in the form of an infinite series, by 
urying out a suitable expansion of the integrand in (5.3), and from this 
s derived the result 


uw" } . >, J(ur) du on 
= sin sé | exp(—vu't) . (5.5) 
WW 6 ow : u 
n 0 
0 
here s = (2n+1)p. (See also Jaeger (6), and Carslaw and Jaeger (7), 
338, ) 
cae, r Pee 
Writing > (5.6) 
2(vt)? 


s last expression may be written in the form 
2 ~ (4s) 1/1 ° 
l sin s§ ——*—_ »* , F, (48; s+-1; — y?) 
i | Loa I'(s+1) 
} 0 


ae (3s) —" P an 
ca. sin 80 M2 1 7 exPl n*) Fi ($8s+1;8+1; y?), (5.7) 


irom equation (2) of section 13.3 in Watson (9). 
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In addition to expression (5.7), which gives the velocity in the form of ap 
infinite series, it is useful to consider expressions which lead to closed forms 








} 
| i} 
| t 
|; 

L M 

Fic. 1. The contour A’ in the «’-plane (Carslaw). 
A’ consists of the curved parts*Aj and A}. 

for w, and in particular for the skin friction 7. A simple change of variable 
in the integration with respect to 6’ in (5.2) leads to the result 


6 
w* — W [ f(n,¢) dd, 
6—B 
where F 
l ' > 79 al , , - 
f(y. 4) - - | exp{—(?+ ’?) F(E,d)n’ dn’, (0.8 
278 . 
0 
E 27’, 
" exp(E cos a’) da’ - 
and F(t.) — | CXP(e cosa ) aa (5.9 
” (E,9) | sin p(d+a’) 


, ¥ 
Finally, since w* = 0 when 6 = 0, the last result for w* may be written 
a 
w* — 2W | f(n, 4) dd. (5.10 


0 
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6. The expression for the velocity in closed form in the case 

B = $n (p =) 

It is of interest to consider an expression in closed form for the velocity 
in at least one case; the case 8 = 37 is of some importance (flow along 
the exterior of a rectangular corner), and it lies midway between those of 
Ravleigh (3) and of Howarth (4). The choice of a special case has been made 
to condense the discussion as far as possible, but it will be clear from what 
follows that the same process would be applicable for any other value of £. 

The contour A’ is now replaced by the two infinite straight lines LL’, 
W'M of Fig. 1, together with small circles surrounding those poles of the 
integrand of (5.9) which lie in the interval —az < a’ < z. In view of the 
symmetry of the velocity distribution about the plane 6 = jz, the dis- 

ussion can be restricted to the range 0 < ¢ < 37. 

The distribution of the poles leads at once to the consideration of the 

following ranges of values of ¢: 
(i) 0<¢ < hdr, 

(ii) @ } T; 

(111) ln<op<n. 
When ¢ = 42, there is a pole at a’ = z, but this we may exclude by in- 
denting the line M’M at this point. Hereafter in this section the classifica- 
tion (i), (ii), (iii) will be used to denote expressions appropriate to the above 
ranges. 

The appropriate expressions for F(£,¢), in each case, are therefore 

; 
(i) F(é,¢) = 3riexp(€ cos d)—2i | exp(—& cosh x)h(x, d) dx 
0 
li) F(&,d¢) 37i — 3771 exp(—&) 
21 | exp(—E€cosha)h(x, 32) dx 7, (6.1) 
0 


(ii) F(E,d) = 3mtfexp(é cos d)—exp(—& sin d)}— 


—2i [ exp(—€ cosh x)h(x, dh) dx 











0 
where 
, _ sin 3(¢-+-77) ; sin $(7—¢) 
h(x, d) cosh $x) ———.— — se7 tf SS re 
sin? 2(¢-+-7)+sinh? 2a © sin? 2(7—d)+ sinh? 2a 
3 3 3 3 
(6.2) 
and h(x, 7) is defined, for x > 0, as 


Osh 2y ain le 
hn, Je) =< eee 
dl sin217--sinh2 227 
sin® 47+-sinh? $a 
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If we now substitute these values for F(£, 4) in (5.8), we obtain a repeated 
infinite integral. Bromwich’s conditions (Bromwich (10), p. 503) permit 
the inversion of the order of integration, and the integration with respect 
to 7’ can then be carried out to yield, ultimately, 


x. 


° . ] l 
i) fn.) = 554.9) +55 [ x(n. evita, 4) de 
“* ] 1 x. 
(i) Sm de) = Fyn —ertii +3 [ ax. 2yh(a. ber) dx 
- * (6.3 
ie 
(iii) f(y.¢) = 9.4 n¥(9, 6) —55 n(n, 7+) 
l x 
a | x(, x)h(x, b) dx 
0 
where w(7, ) cos d exp( 72 sin2)f{I Lerf(7 cos 4)} . 
x(7), #) cosh x exp(7* sinh*x){1—erf(7 cosh x)} | 
2¢ 
- erf(y) = = | exp(—z*) dz. 
0 


Before the final substitution is made in (5.10), it is necessary to consider 
the infinite integral ; 
G(n, d) x(n, x)h(a, hb) dx. (6.5 
0 
The discussion of G presents no point of special difficulty, and it is sufficient 
here to state the results. Taking » > 0, the integrand in (6.5) is a con 
tinuous function of ¢ and 2, if ¢ is restricted to lie in either of the ranges (i 
and (iii). Under these conditions, G(y,¢) is uniformly convergent in ¢ in 
these ranges, and is therefore a continuous function of 4. Furthermore, it 
can be shown that G@ is bounded throughout the region 
ln—e <¢: Ilm@te (€ > 0). 

The properties of G permit the inversion of the order of integration in the 
repeated integral which arises when the expressions in (6.3) are substituted 
in (5.10). The integration with respect to ¢ can be carried out to give, 
when 0 < 6 - 


Lar 


D) 


w* 


5 erf(7 sin 6){1-+-erf(7 cos 6) + 


l . , ‘ i 
t= n sin d exp( 1? cos*d)erf(n sin $) dp-+— nx(n, x)k(x, 8) dr, 


. 


0 0 (6.6) 
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where 





k(a, 0) = log ! —— 
7 \ 


‘cosh 3a— cos 3(6-+-7)Keosh 3a+-cos 3(7- m1) (6.7) 
_ vealed 2. hoses ned? 4 5.7 
Sa 2y 2 ar\\Seosh 27 HS 2 

‘cosh 2a-+-cos §(0-4 7) cosh 3° — COS §(77- 0)} 


The corresponding result for w* in the region $7 < 6 < zm is obtained 


y separating the integral in (5.10) into three parts, namely, 


i 7 @ lorte 
7 2 [ fdd+2 [ fdd+2 [ fad, 
where 0 € @—4n. 


The last of these integrals is of order e, and clearly integrals of the type 


sir —€ 
us dd 
0 
iv be replaced by us dd Ole). 
Thus 
* r j rl ’ 1 
7 NY, @) dd | ruh( 1); 9 7T T dh) dd 
. 


nx(n, x)k(x, 0) da 


nx(n, x){k(a, Jr—e)—k(ax, +-€)} dx+-O(e). 


It is a simple matter to show that the last integral here is of order ¢, while 
the first two integrals combine into a simple form. Finally, therefore, since 


«may be chosen arbitrarily small, the expression for w* in the range 


T H 7T 18 
_ } erf(7 sin 6){1+-erf(7 cos 6)|— 4 erf(7 cos @)-4 
3{1—exp(— 7?)} — x(n, «)k(x, 0) dx. (6.8) 
0 


the usefulness of the expressions (6.6) and (6.8) is restricted in that the 
ninite integral contains the two parameters 7 and 6. However, for 
suthciently large values of 7, the asymptotic expansion 


erf(n) ~ 1 exp | Oh : l Ae 3 7 15 + -| 


27?  4n* 878 


aie 
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leads to the approximations 


toh 


27! . 


exp(—7”) (0 <6 - 
0 sexp(—7?) (40 < 6 - 


| - eC bn) ) 
<= | ax(n,x)k(a, 0) dx ~ |x i sd 
7), 


Furthermore, there is an approximation for the first integral occurring 
; in (6.6). Writing 
4 
l . ea , 
A(n, 0) = = [ n sin d exp(— 7? cos*d)erf(7 sin 4) dd, (6.10 
0 
it can be shown that 


] 
0 < }{erf(n)—erf(y cos 0)!|— A(n, 0) < = nexp(— y?){1—cos 6}. 


(6.11) 
Summing up briefly, therefore, the expressions for w are 
Fi = 1— Serf(7sin @){1+ erf(7 cos 4)! 
0 | 
] . | 
— = [ 7 sin ¢ exp(— 7” cos*¢)erf(n sin d) dd— | 
a (6.12) 
0 e 
5 | 
(valid in the range 0 < 6 < 4m) 
7 = 1— Serf(nsin 6){1+ erf(7 cos 6)}+ 3 erf(7 cos @)- | 
i. | 6.13) 
${1—exp(— n?)} — = } x(n, x)k(a, 0) dx aac 


(valid in the range }7 < 6 < =) 


These expressions are sufficient to describe the motion throughout the 
region occupied by the fluid, in view of the symmetry of the motion about 
the plane 6 = 37. For large values of 7 the approximations (6.9) and (6.10) 
can be used for the integrals in (6.12) and (6.13), and the expressions there- 
fore provide general information in the region where the convergence of the 
solution in series would be slow. Fig. 2 shows curves of constant velocity 
for the case 8 = 3x, drawn from tables given by Jaeger (6) for the corte- 
sponding heat conduction problem. 
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O 0-5 10 X 15 

‘1g. 2. Curves of constant velocity near the corner in the case f Sor (p z). 

Fic. 2. ¢ f | tl Sr ( 5) 
In the diagram X 7 cos @ and 7 r/2(vt)3. 


7. The skin friction 
Denoting the skin friction by 7, and considering the friction on the plate 
0, then (for ¢ 0), 


[ow B ow 7 
T ~ Daf )bar ; ‘ 
Saal, areas), 0 = 


0 
The expression for 7 as an infinite series is obtainable at once from (5.7). 
Term-by-term differentiation with respect to @ is permissible, since the 
series is uniformly convergent in 6, and so is the derived series (for any 


range of V alues of G0). Thus 


; scar “ tT = 7° (38; 8+1; —7?), (7.2) 
inequivalent »>%rm, 
= > exp s\n a, (7.3) 
B(vt)'n 47°} u 
where 8 (2n-+-1)p (2n-+-1)z/B. 


From (7.2) it appears at once that 7 takes the form 7”-» as 7 > 0; that is, 


is not bounded near 7» = 0 when 8 > z, though (as will appear later) the 
irictional force on any finite part of the plate is bounded (for t > 0). 

It will be necessary later to consider expressions for 7 for particular 
values of p. We can indicate these cases briefly by introducing the suffix p, 
when necessary, on 7 and on any other quantity we may wish to consider. 
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Thus, for example, 7, will refer to Rayleigh’s value of + for the flat plat 
problem; that is, uW 
fs (avt)* 


We consider next expressions in closed form for 7. From the result (5.10). 


6 
w 


7 = 1-2 | fd) ag, 


fOw - 
so that (- 2Wf(m, 9), 

og 6-0 
since f is continuous in ¢ at d = 0, provided that p is not an integer. Hene 

ce . 

F— f(y, 0). (74 
(vt)! » 
The evaluation of f(n,0) proceeds exactly as before. The contour 4A’ js 
replaced by the straight lines LL’, M’M, and small circles surrounding 


those poles of the integrand of (5.9) which lie in the range —z < a’ <7 
Since the case where p is an integer has been excluded (for the moment 
let n’ be the greatest integer such that n’ < p (so that n’ = 0, if p<] 
Then the poles lying in the given interval are the points 


, 


x sB, where s 0, +1, +2 Ln’, 

As in section 6, we derive a term which is a repeated infinite integral in the 
expression for f(7, 0), and an appeal to Bromwich’s conditions again permits 
the inversion of the order of integration. The integration with respect to: 
may be carried out to give, after some reduction, 








. = . sin px [ y(n, x)cosh pa da 
r= 4,[Hl+erf(q)}+ ¥ (—1)6(n, 98) 4+ rk | Xiu Zoos pee 
- s=1 B J sin®p7+sinh*pxr 
0 
(7.9 
This is valid for n’ > 1. If n’ 0 (p < 1) the expression is 
, sin pz [ y(n, x2)cosh pa dx che’ 
T = 7,|${1+ erf(y)} : ZA ~ ; = (4.0 
= B sin*p7+-sinh* px 
0 . 


These two expressions deal with all cases, except where p is an integer. 
We could proceed directly from (5.10) to discuss this case, but it has been 
dealt with in a previous paper (Sowerby (5)), and all that seems necessar\ 
here is to state the result. 


That is, for p = n, where vn is an integer, 


7 4{1 +-erf(y)}+ 3(—1)"+4{1—erf()}4 > ( 1)*(7, 88)}. (7.7 
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if we consider further the case where p < 1, we observe that the infinite 
steoral in (7.6) converges only if » > 0 (as expected). However, the 
iivergent part may be separated from the integral as follows: 


cosh x cosh px ss 
Since — —~— = 2cosh(1—p)zx+ U, 
sin“ pi sinh*pa 





, (3—4sin2pz)cosh(1— p)a—cosh(3p— 1)x 
- 2(sin*p7-+sinh?px) 


y(n. x)eosh pa da a . 
| r 2 | cosh(1—p)x exp(n? sinh?x){1—erf(7 cosh x)} dx 
7-+-sinh*pa 


sin-p 


( U exp(7* sinh?x){1—erf(7 cosh x)} dz. 


0 
Yow | Oi where A > 0 (since p 1), so that this last integral 
nverges for all 7 0), Write 3 | P) m(Q<<m |). and consider the 
rst integral 
2 | cosh 2ma exp(7* sinh*a erf(7 cosh x)} dx 
cosh 2mx exp(7*? sinh*x) dx | exp(—2*) dz 
0 neosha 
{ . F 
exp| > cosh 2ma dx | exp(— y?—2ny cosh x) dy 
0 
t 


exp(—7*) | exp( y?) Ky (2ny) dy 


“mn 


sec M7 exp(— +s 7") K,,(3?). 
It is not necessary here to justify the change in order of integration, since 
his result could have been established by proceeding directly from the 
xpression for f(m.0). The last result is deduced from Hankel’s integral 
Watson 9 pp 393-4). 
Finally, therefore, the expression (7.6) may be written 


») 


T 7; l erl(7); — sin” 7 eXp( by?) K,, (377) + 
U exp(7* sinh*x){1 -erf(7 cosh x)} dx (7.8) 


(p l, and y(1 Pp) m). 
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Finally, results (7.5)-(7.7) may be used to derive asymptotic expansions 
for r. For example, selecting the case where p < 1, and using the result 


exp(— 7?) ] 3 15 
x(7,@) © see ee +, 


yn 27? cosh2x § 4n*cosh42 8n* cosh®x 
] 7) 1) q 


we have, from (7.6), 











1 alCs , & , & ) i 
Tw 1+ — exp(— 7°) —=+4+-44+---}], (7.9 
8 ie a 
where Cy = tas, C, = —Zas, C7 = Ba, 
2sin pr [ cosh px dx - 
and a, = 1—-—~— | — — — . (7.10 
B J (sin’pw+sinh?p2)(cosh x)"—! 


0 
(There is no term of the form exp(— 7?)/7 in (7.6), since «, = 0 for p <1) 
The coefficients «, can be evaluated whenever p is a rational number, and 
the two cases p = } (Howarth’s case of the semi-infinite plate) and p = ? 
(8 = $7) will be considered. In the former case the «,, are conveniently 
evaluated by contour integration; in the latter, the integrations are more 
tedious. Evaluating the first few terms in each case, we obtain 








v . 1 x 2) l 4 § ; 2265 | 
. ~ Se = 2X _ ad — pr ae 2 lee 
Ty Ty = = ” 873 128° | 409677 | | (7.11 
. ts 
7 oo af 4 64 | 
le ae a att] 





The expression for 7, differs from the result (80) given by Howarth (4). 
Professor Howarth agrees that there is an error here (private communica- 
tion), and the above expression is in agreement with his result (79). 


8. The frictional force parallel to Oz 
The frictional force on a strip of plate of unit length and width d, extend- 
ing from the corner, is a" 


F [ 7 dr. 


0 


Thus, if F, is the force in the flat plate problem, 


Wb 
F=f. 
1 (avt)! 
and b b’ 
F—F, = | (r—7,) dr = 2(vt)! | (r—7,) dy 
0 0 
T —2(vt)! [ (r—7,) dy, (8.1) 


b 
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he b 

where , ’ 

2(vt)! 
nd T 2(vt)? (r T,) dn. (8.2) 

0 
Ty ‘nite integrals here c roe. whatever the value of ior 
The infinite integra ere converge, whatever the value ol p (p = +), since 
- is of order 7”~) near 7 = 0, while t—7, is at most of order exp(—A7?) 


for large values of » (A depends on p, and is > 0). It follows also that if b’ 
ssufficiently large the last term in (8.1) will be very small, in which case 7’ 
will represent closely the excess of F over F,. 

We proceed therefore to evaluate 7’. This could be done by substituting 
the expressions for 7 ((7.5)—(7.7)) in (8.2), and inverting the order of integra- 
tion in the resulting repeated integral to give, ultimately, expressions for 7' 

closed form. Rather surprisingly, however, the result appears in its 

eatest form if we take instead the result (7.3), where 7 is expressed as an 


nfinite series 





Thus 
2 ‘(2m Ff *\. 
T= Wi > 8 | exp(- ve) du) dy 
J [sna , 4y7] u 7) 
uW | {B | | dn, (8.3) 
vhere s = (2n+-1)p, and 
B . > exp| pe 
a=o. 3 1 - 
lee u? 
aD, | EXP gaa} ea + Spanley} 
Pl n=0° £ 
~ > exp(—4yTy-vb0?) + yoiv(h0)}, 


0 


> ‘ 


Watson (9), equation (1) of section 3.2 and equation (5) of section 13.3. 
In order to perform the summation, we express the Bessel functions J in 
terms of contour integrals (Watson (9), equation (3) of section 6.22). We 


then obtain, placing z Ly? 


exp(z cosh w—4sw)cosh }w dw. 









este ws 8s ee * 
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Interchanging the summation and integration, we have, on summing the 
series of exponentials, 


B l “exp{2( (cosh w—1)} cosh }w dw 
2728 . sinh(4pw) 
provided |e”) < l,i.e.re(w) > 0, so that the contour must be to the righ; 


of the imaginary axis in the w-plane. 
The reversal is legitimate, as it is easily shown that 
| ‘exp(z cosh w) || cosh }w > exp(— }sw) dw 
converges along the contour. 

Since the only poles of the integrand are at the points w = 0, w + 2nig 
where v is an integer, we may deform the contour up to the imaginary axis 
indented by semicircles of radius 5 round and to the right of the poles. W, 
shall take for contour the straight lines joining the points «—iz, —iz 

im, ©-+im with the imaginary axis indented as explained, and then let 
o—> 0. 

The integrals along the straight lines on the imaginary axis cancel one 


another in pairs; the integral round the semicircle at the pole w = 0 is 


easily seen to be 1/7!; the sum of the integrals round the poles w = +2unif 
amounts to m o 
™ — (—1)" cos nB exp{z(cos 2nB—1)}, (8.4 
— “7 
n=] 
where m is the integer defined by the inequality 7/28—1 < m < 728. 
This still applies if 7/28 is an integer, since the points corresponding to this 
— of n, namely w Lim, are not then poles. Further, if 7/28 < 1, ie 
B 7, this term disappears. 
We : are left with the integrals along the lines parallel to the real axis 
For these we write w im+x and w = im+2 and obtain 
2 i  % » >) > ee, > 
Bb I+ > 2(—1)" cos nf exp(— 2zsin?nB) 
7 n=] 


x 


tcos$pz [ sinhaxsinh pa exp(— 22 cosh*x) dx 


B cosh 2pa—cos pr 


Hence, from (8.3), we have 


" 4 W ~ 
T — 7b (5 1)" cos nB exp(— 2z sin?nf)-4 
W* n 


x 

, 2cos$pz [f sinhasinh px exp(— 2z cosh*2) dic| dn. (85 
a - = om _ —<$—$_— 7). Out 
cosh 2pa—cos px 
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Hence, reversing the order of integration in the last integral, and using the 


| 
Suit 


l 


“y 


} exp( y”) dy 
e obtain finally 
’ » 1 “é” ot ° . 
' 2 COS 5 p7r tanh xsinh px da : 
7 2uV S { 1)” eot nB - at, inant P (8.6) 
=! 3 | cosh 2px—cos pz 
. 0 


? 


the reversal is legitimate provided the integral in (8.6) exists, its integrand 


ing always positive. 
un 
This result applies in the general case 8, bearing in mind that the term > 
1 


inishes if 6 lr, as explained previously. 
The integral in (8.6) can be evaluated in closed form for B ar /1, where 
und / are integers such that 0 B <2n. If B a/n, where n is an ' 
ntegel 2, a simpler form of the result is available, starting from the y 
xpression (7.7) for r. In the resulting integral for 7’, differentiation of each 
term with respect to 8, and an integration by parts with respect to 7, vields : 
he result j 
ull nal : 
7 1+(—] ‘1.9 > (—1)8(~— sB)cot sB}. (8.7) ’ 
a al b 
lable I gives the values of 7’ for selected values of 8. The value for B = 27 ‘ 
isalready been given by Howarth (4), and the table extends other results 
en by Sowerby (5), with corrections in the last column to two mistakes 
nhis A; and A,). We note at once the rather surprising result that, if fluid 
ows along both sides of a right-angled corner, the total drg is less than y 
he corresponding drag over both sides of a flat plate. 
9, Discussion of the results 
\ 


Chere are two ‘boundary-layer thicknesses’ involved in the wedge prob- 
m. The first is the usual one in connexion with the flow of fluid over a flat 
ite, which we may call the ‘main’ thickness; the second is the width of 
he region in which the effect of the corner is significant. We may define 
he width of this region to be the distance from the corner of the point at 


vhich + first assumes the value 7, to within 1 per cent. 


Firstly, the width of this region is of the order of the main boundary- 


yer thickness. For example, if we examine the case 8 = 37 and write 
‘= ncos0, Y 7 sin @, it is evident at once from result (6.12) and the 
pproximations (6.9) and (6.11) that when X > 2 the flow is not dis- 


tinguishable from that over a flat plate (in fact, the statement is true if 


A l-1). The thickness of the main boundary layer is Y 2. 
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TABLE | 








B TW T pW 
Qa 4 0-5 
8v3 2 

Sar a 0-3010 

5 2.06 oO 

7 0 0 
l6v3 2 

ar a ; 0-2386 
ai o7T 


har 0-6366 
7 
9 
x = 1-155 
iw € 
V3 
9 | 
{7 —1 | —1-637 
7 | 
Lr 2 cot $7+ 2 cot 3a 2-103 
2 10v3 = 
$7 ; 2-561 
7 9 | 
Lar 2 cot 47+ 2 cot ?7—2cotix 3-015 





Secondly, the width of the region changes with the angle of the wedge. 
For 8 = 27, Howarth found the 1 per cent. point for |+—7,|/7, to be where 
7» = 1; for 8 = 37 and }z the corresponding results are 7 = 1-1 and 1:85 
respectively. Other values have been given by Sowerby (5) for smaller 
values of 8, and it appears that the width of the above region increases as £ 
decreases from the value 27. 

Further, considering the whole boundary layer in the general case 8, the 
thickness decreases as we approach the corner if 8 > 7, and increases if 
B < 7; these effects can be seen respectively in Fig. 2 and in Carslaw and 
Jaeger’s figure for the rectangular corner ((7), p. 151). 

The effect on the skin friction 7 can be seen from result (7.2), and the 
asymptotic expression (7.9) for 8 > 7 (a similar expression is available for 
B <7). Thus for 8B < z, 7 decreases steadily from 7, to zero at the cornet 
while for 8 > 7, 7 increases steadily from 7, and is not bounded at the 
corner itself. 

We consider finally the frictional force per unit length when a rectangulat 
pipe of sides a, b (a > b) is set into motion in the same manner as the 
original boundary, the liquid now filling the interior of the pipe and also 
the whole space outside it. We shall require some restriction on } to ensure 
that the corners do not interact with each other. The force on the interio! 
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fsuch a pipe has been determined previously (Sowerby (5)) and the 


»striction on 6 which is sufficient is that 
b 8(vt)?. 


\s we have seen, this is also adequate for conditions on the exterior of the 
pe (where 8 = 37), and it is a simple matter to show from (7.11) that the 
st term in (8.1) is certainly negligible in comparison with 7’ when 6 

aatisfies the above condition. Using the appropriate values of 7' from 
lable I, we obtain for the frictional force F per unit length (taking into 
count both the interior and exterior faces of the pipe) 





, if 
P = 4.W|—_— al, (9.1) 
(arvt)* ; 
here A : (T,-1-T.) 0-6712. 
uW : 


Ve may now take the modified form of Rayleigh’s hypothesis, as explained 
previously, to estimate the drag on a pipe of rectangular cross-section with 
sides a, b (a > 5b) and of finite length /, the undisturbed velocity of the 


steady incident stream being W in a direction parallel to the length of the 


pe. Thus if z is distance measured downstream from the leading edge of 


the pipe, we replace ¢ in the above results by y?z/W. The restriction now 


mposed on 6 is (writing z l) 
. (vl \4 
b 13-6 aie (9,2 
WW 
nd the force F per unit length is 
; 1+b/W\h 
F = 4u0 | ( )°. al. (9.3) 
y \avz 
Hence the drag coefficient C,, for such a pipe would be, taking into account 
| faces. 
l ee 2 2-685v 
Cy F dz AE. es (9.4) 
pV “Hla b) B R (a + b)W 


where a 1-32824 and R Wliv. 


These results, of course, are subject to the limitation (9.2); that is, they 


re valid only for a suitably short length of pipe. 


Une of the authors (L.S.) wishes to express his gratitude to Professor 


hosenhead and Professor Howarth for their kind advice and encouragement 
luring the course of his work. 
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Note added 17 November 1952 


After this paper was accepted for publication the authors were 


informed of the recent publication of a short note by Hasimoto (11 
Hasimoto obtains a table of values of 7'/uW for eight values of 8 in th 
range 47 <fB < 2z, and the corresponding results of Table 1 are in 
agreement with his values. 
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(HE TWO-DIMENSIONAL IRROTATIONAL FLOW OF A 
COMPRESSIBLE FLUID AROUND A CORNER 
By H. J. DAVIES (The University, Southampton) 
Received 31 July 1951] 
SUMMARY 
\n exact solution is derived for the two-dimensional, steady, irrotational motion 
mpressible fluid around a corner. This solution provides another example in 


h there exist stream-lines in the physical plane which are free from singularities. 


so provides an example in which Lighthill’s definition of a solution of a hyper- 


tric differential equation is used in the case where the value y 1 is a simple 
f the function F(a, 8; y; 7) regarded as a function of y. The hypergeometric 
ns required in this solution are tabulated. 


|. Introduction 
Usinc the basic assumptions of continuity, irrotationality, and that the 
w is isentropic, Chaplygin (1) obtained the following fundamental 


ations of the hodograph method: 


(1.1) 


. which we obtain 


| 2,/ / 
Ar | T | T 37} | {(1—+)— 2Br 
te CT ' ’ C6? 


) Cc =u 


0: (1.2) 


re 7, 0 are the hodograph coordinates used by Temple and Yarwood (2), 


q? Geax» Where q denotes the fluid speed, ¢,,,, the maximum speed 


the motion } Li(y 


), and y the adiabatic index. 
The solution of equation (1.2) adopted here is that suggested by Chaplygin 
1), i.e bi cid 

ws Ar”? F(z7)sin mé. 


\ ») 


ibstituting in equation (1 it is found that the function F (7) satisfies a 


pergeometric differential equation, the two independent solutions of 


which are 


F(a,b;ce;7r) and 7!-¢*F(1+a—c, 1+6—c;2—ce;7r), 


} ab sm(m |-1)B, c m+1, B 2-5. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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It is seen that when m 2 one of the solutions of the differential equation 
satisfied by the stream-function is 


sb = ArF (2-5, —3; 3;7)sin 20, 


where F(2-5, —3; 3;7) ay(32— 807-+-707?— 217%) is tabulated in Table [+ 
This solution, which represents the flow of a compressible fluid in a corne 
was investigated by Williams (3). The solution is characterized by , 
singular point in the hodograph plane, which is an extremum of 4%, trans 
forming to a ‘point stream-line’ in the physical plane. The value of A js 
chosen such that ¢,,,. = 1; we find A 1-034722. The maximum speed of 
flow occurs on the ‘point stream-line’. This speed, q, is equal to 1-367q, 
where q, is the critical speed, giving a maximum Mach number of 2-262. 


2. Flow around a corner 

The second of the two independent solutions is more interesting in that 
some of the stream-lines in the physical plane are free from any singularities 
Following Temple and Yarwood (2) it can be shown that this solution 
represents the flow of a compressible fluid around a corner (see Fig. 1). 

Another point of interest in this solution is the fact that it provides an 
example in which the value (2—c) 1 is a simple pole of the function 
F(1+a—c, 1+-b—c; 2—c:;r) regarded as a function of (2—c). In this case 
Lighthill (4) defines the solution of the hypergeometric equation as 


F*(0-5, —5; —1;7) im] (03+ 28, 5— $30; —1—8;7) 
50 


o— 


1+ 2-57—7-5F (2-5, —3:3;7)log r— 23-7573 + 36-6015625074 


>» 


5+r)C(r—5y)(r lyr, ( 


- —. [*(0- 
17-433593757°+ 90 S dint. BA ts : 
— P(2-5){P(r+ 1) 


For values of 7 greater than 0-5 the above expression for the expansion 0/ 
the function converges very slowly. A more convenient form in the region 


0:5 <— + < 1 ean be obtained as follows. In the form taken for the stream- 
function 


us 


7-1 F* sin 26, 
the function F*(0-5, —5;—1; 7) satisfies the hypergeometric differential 
equation 2K dF 

z(1 r)- - + fe—(a+b-+ 1)z}- —abF 0. (2.3 
dr* dr 


+ Tables I-IV will be found at the end of the paper. 
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where a = 0-5, b 5, ¢ 1. Substituting 7 = 1—+t in this equation 


ve obtain the hypergeometric equation 
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Fic. 1. Stream-lines in the physicai plane. 


the two independent solutions of which are 


F(0-5 5;—2-5;l—r) and (1—r)*5F(—1-5, 4; 4-5; 1 


which are valid when 0 < + < 2. Thus the region in which t 


solutions are valid has an area in common with the region —1 - 
which F*(0-5, —5: —1:7) is valid. In this common area the three 


must be connected by a linear relation, say 


F*(0-5, —5; —1;7 A F(0-5, —5; —2-5; 1—7)4 


B(1—r)*>F(—1-5, 4; 4:5, 


the other independent solution of equation (2.3) is given by 


72 F (2-5, —3:3:7). 


Hence the solution F(0-5, —5; —2-5; 1—r) is linearly connected, within the 
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common area, with the two independent solutions F'*(0-5, —5; 
7? F (2-5, —3; 3:7) giving 

F(0-5, —5; —2-5;1—r) = A’ F*(0-5, —5; —1;7)4+ B’r?2F (2-5, —3: 3:7), 
When 7 = 0 


A’ F'(0-5, —5; —2-5; 1) - 
1'(2-5)I 





r(—2-5)P(2) 
¥..3) 


Hence we may write 
F*(0-5, —5; —1;7) = A"r2F (2-5, —3;3:7)- 
+ B(1—r)*5F(—1-5, 4; 4-5; 1—7). (2.4 
When 7 = 0 we get | BF(—1-5, 4; 4-5; 1) so that 
l 1'(0-5)1(6) 128 


———— - 
F(—1-5, 4; 4-5; 1) 1(4-5)P(2) 7 


When 


7 l we get 
1" ll a 5) “ 1; 1) ai 
F (2-5, —3; 3; 1) 
where F(2-5, —3:3:1 1(3)1 (3-5) ; | 


P(6)P(0-5) 32° 
From the definition of F*(0-5, —5; —1;7) given by equation (2.1) we obtain 
Paro, —G: —]< 1) 














lim} F'(0-5-+-335, —5—338; —1—86; 1) 
5—0 
1S (2-4-3) .. Pere i 
ane ) F (25 |-355, —3—138; 3-+-5; | 
4 ra) a ae 
R r'(—_1—8 Lh:(2-+3 1(3-+-8 
I'(3-5)lim | —— A. aoe ~ a ~ ——— Bh in 4 rT4 
so] ['(— 1-5 — 335) (44-385) 4 & LP(0-5—383)P(6+330) 
Substituting 
l aon (0-5) |? 
I'(—1—s) = ———_——,, and [I(—1-5—8§) BEX i. 
61(2+8) = ['(2-5-+-336) 


and using the relationship 


x 


A. tied Se ea wn ee | 
I'(z) i 2 — n(n-+-z) 


where y is Euler’s constant, on proceeding to the limit we obtain 
r*(0:5, —6: —1; 1) 43(3—2 log 2). 
Substituting this value in (2.5) we obtain A” 7:5(1-6137056), so that 
from (2.4), we have 
F*(0-5, —5; —1;7) 1-6137056 x 7-572.F'(2-5, —3; 3;7)-4 


4228(1 7) F(— 1 
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Thus the function #'*(0-5, —5; —1;7) is evaluated for values 0(0-1)0-5 of 7 
rom equation (2.2), and for values 0-5(0-1)1-0 of t from equation (2.6)—see 
lable II. 


3, The critical stream-line and maximum velocity 
Craggs (5) has shown that the curve 
O(x,Y) 
o(T, @) 


0 
represents, in the physical plane, the cusp locus of the stream-lines. This 
urve transforms into a curve in the hodograph plane whose equation is 


O(dbs, hd) 
o(7, @) 


Substituting for % and ¢, this equation takes the form 


) 1/6 
tan 20 r |= z} (3.1) 
~tP’ a \l—r 
vhere P= ge" 
nd r dP /dr. 


[his curve divides the hodograph plane into regions, each region, when 
taken separately, transforming into a physically admissible flow pattern 
nthe physical plane, the transformation being one-to-one. 

Craggs (5) has shown that if an extreme value of % does not exist, i.e. 
P’ + (0) (as in this example), then there will exist in the hodograph plane 
stream-lines which do not cut the boundary curve. These stream-lines will 
be bounded by that stream-line which touches the boundary curve. Such 
1 stream-line will be referred to as a ‘critical’ stream-line. 

Let the equation of this stream-line be 

ub AP sin 20; (3.2) 
eliminating @ between equations (3.1) and (3.2) we obtain an equation for 7 
»f,2 P4 67 ] 
uJ (07 ) 9 « 
; i —; (3.3) 

A? 7? P’?(1—7)+ P?(67r—1) 

which gives the points of intersection of the stream-line with the boundary 
curve. As (3.2) is the equation of the critical stream-line this value of 7 will 
be a double root of (3.3) and hence will be a single root of the total derivative, 

e.a single root of the equation 
} (7) 2(67 1)2P2 1573 PP’ 2(67 1)(7 1)r?P’2 Q). 

By computing values of f(r) together with the first, second, third, and fourth 
differences and using Bessel’s interpolation formula, it is found that f(r) = 0 
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when + = 0-22524. Hence the critical stream-function is obtained by sub 
stituting this value of 7 in (3.3), giving 


b/A 4-6897. 


Substituting this value of 4/A and @ = 45° in (3.2), the value of P corte. 
sponding to the maximum value of 7 associated with those stream-lines 
which are free from any singularities in the physical plane is obtained, viz, 


P = 4-6897, 
which gives, by graphical interpolation, 


Tmax 0-355, 
4. Equations of transformation 


In the general case when 
ub Ar”? F(7)sin m6 


Ringleb (6) has obtained the transformation equations. With reference to 
these, taking m = 2, the required transformation equations are 


fas 2 7~*(1—7)-*5[ 3(7P’ — P)cos 30+-(7P’ 4 P)cos 6], 


imaty 7-*(1—7) 2 (7 P’ — P)sin 30— (7 P’ + P)sin 6]. 
5. General discussion 

The flow of a compressible fluid around a corner may be described as 
follows: for the entire flow, the velocity is zero at the entrance and at the 
exit at infinity. For a closer examination of the flow we divide it into two 
regions which, in the hodograph plane, correspond to: (1) stream-lines that 
do not cut the boundary curve, and (2) those that do. These two regions 
are separated by the critical stream-line y:/A = 4-6897 which touches the 
boundary curve (see Fig. 2). 

Type (1) transforms into stream-lines free from any singularities in the 
physical plane. Along some of these stream-lines the fluid velocity rises 
from subsonic to supersonic values and returns again to subsonic values. 
[t is seen in Fig. 1 that such stream-lines form paths which converge from 


infinity, reaching a minimum section of the path in the neighbourhood of 


the line where the speed is equal to the local speed of sound; the stream-lines 
then diverge to the line of symmetry. The maximum speed occurs on 
the critical stream-line, /A = 4-6897, and is given by + = 0-355, i.e. the 


maximum speed, qg, is equal to 1-46¢, 1/(1+ 28)! .qmax is the 


a) 


where q, 
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ritical speed, giving a maximum Mach number of 2-75. Along the critical 
¢ream-line two irregularities appear at points at which t = 0-225. These 
~egularities of the critical stream-line are also singularities (cusps) of the 


usp locus (see Craggs (35)). 


90° 80° 


70° yw =Ar'F’ sin 20 
60  F—Table II 
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Fic. 2. Stream-lines in the hodograph plane. 








Beyond this critical stream-line we have the second flow-region. On every 
stream-line in this region the maximum speed is greater than that given by 
r= 0355. Also, four irregularities appear in the form of cusps at points 
corresponding to intersections of the stream-line and the boundary curve 
inthe hodograph plane. As is seen in Fig. 1, this region presents a physically 
impossible pattern of intersecting stream-lines. However, if the corre- 
sponding region in the hodograph plane is regarded as divided into separate 
sub-regions by the boundary curve—these sub-regions being indicated by 
continuous lines, sectional lines, dimensional lines, and again sectional lines 
in Fig. 2—then each sub-region, when considered separately, presents a 
one-to-one transformation to the physical plane. Hence each of these 
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regions presents a family of stream-lines which is theoretically acceptable 
Of these families, the one represented by continuous lines, being of subson; 
type, i.e. acceleration and hence negative pressure gradient occurring in thy 
direction of convergence of the stream-lines, has a correspondence with tly 
incompressible case. 
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THE SOLUTION OF DIFFERENTIAL EQUATIONS IN 
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THREE DIMENSIONS 
[I. POTENTIAL FLOW ROUND AEROFOILS 


By D. N. DE G. ALLEN and §S. C. R. DENNIS 
(Imperial College, London) 
Received 18 September 1951] 


SUMMARY 
Part 1 of this series (1) described how the relaxation method could be extended 
two to three dimensions in order to solve boundary value potential problems. 
his paper conside rably more advanced solutions of Laplace’s equation are pre- 
1; they require the introduction into three-dimensional relaxation of graded 
es and of the determination of initially unknown boundaries or interfaces. The 
lems solved in illustration concern the potential flow of an incompressible fluid 


indthree-dimensional aerofoils ; two cases are solved in detail, one for the flow round 
ght wing and one for the flow round a swept-back wing. Results are presented 
form of numerical values of the velocity-potential of the flow at the points of a 
xation lattice, together with certain derived curves which show the variation of 
culation and of the position of the aerodynamic centre along the span. The 
m in three dimensions necessitates a formulation which would not normally 
ide for the corresponding problem in two dimensions ; consequently the method 
lescribed firstly, for simplicity, in relation to such a two-dimensional problem. 
xistence in this case of an exact analytical soiution makes possible a comparison 
relaxational results. 


Tue two-dimensional flow of an incompressible fluid round an aerofoil 
iid normally be solved in terms of the stream-function y; if axes of x and 
re respectively parallel and perpendicular to the direction of undisturbed 
w, the corresponding components of velocity, u and v, at any point in the 
Ww are respectively given by 

CW . Cus 


cy CX 


vhere u satisfies L iplace’s equation 


V2u; C= \ O7us 0. (2) 


ox" oy" 
gether with the conditions that, at large distances from the aerofoil, 
Cu Cus 
» e/. > 0, (3) 


Cy Cx 


ere LU’ is the velocity of the undisturbed flow; also that, on the boundary 


he aerofoil us k, (4) 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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where k is a constant whose actual value in any particular case is determine) 
by some extra single assumption or postulate: this extra condition woul 
usually be that the stream-line % = k leaves the aerofoil tangentially at th 
trailing edge. 

Stated in this form the problem in two dimensions offers no difficulty t 
solution by relaxation methods; it can be found by means of a synthesis 
between two numerical solutions of equation (2), which satisfy the con 
ditions (3) and also (4) for any two different but arbitrarily chosen values 
of k; and when ¢% has been thus determined, the components of velocity an 
given immediately by (1); other quantities, such as the circulation round 
the aerofoil or the pressure distribution round its boundary, can be caleu 
lated directly. 

2. A statement of the problem in terms of % in two dimensions cannot 
be extended to three dimensions since in three-dimensional flow no strean 
function exists, and we are compelled, therefore, to use the velocity-potentia 
¢. A new technique is therefore essential in three dimensions; it will by 
described, for simplicity, in two dimensions, where it is not essential and 
where indeed it is not so suitable as the technique suggested in section 1. 

In terms of ¢ the components of velocity are given by 
Cd . Ch 


Ox oy 


Uu 


and the velocity potential-itself satisfies Laplace’s equation 


V*o = 0. b 
The boundary conditions become 
c . c : 
g is p > 0, 
Cx cy 


at large distances from the aerofoil, together with 
Ch 


CV 


Q, (5 


on the boundary of the aerofoil itself, where @/év denotes differentiator 
along the normal to the boundary at any point. From the point of view 
of arelaxational solution there is little to choose between (6) and (7) in terms 
of ¢ and (2) and (3) in terms of ys; but (8)—a boundary gradient condition 

is less easy to satisfy than (4)—a boundary value condition. A more impo! 
tant feature is that, unlike ¢, d is not a single-valued function; it is instead 
eyclic, for the circulation, K, round any closed contour, C, in the plane 0 


K » ds 
as 


Cc 


the flow is given by 
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rmine where s is arc-length measured round C), and K does not vanish if C 
wou encloses the aerofoil. 
att 3. Fig. 1 illustrates how this last difficulty may be overcome: 
ult / e 
nthe 
1e ( 
>_> 
~ 
| TO 
ca 
canl 
tre r 5 
tent Pp ) 
will 
ial 
ion | y 
x 
G H 
Fic. 1. 
i {Brepresents an aerofoil placed in a fluid stream which, if undisturbed by 
the presence of the aerofoil, would flow with uniform velocity U’ in the 
direction. BD represents that stream-line of the flow which leaves the 
» | wailing edge B of the aerofoil tangentially. We imagine the plane of 
the flow, which is actually cut by the slit AB, to be also cut by the stream- 
tiatir ne BD, and consider the region of integration of the governing equa- 
of vit tion (6) to be singly-connected and bounded by the external boundary 
terms | ABDEFGHDBA,; the line ABD is thus a double part of the boundary,t 
ition ind at all points along it the velocity-potential ¢ has two values, one 
impol elonging to the region above A BD and one belonging to the region below: 
nstea to distinguish them we denote ¢-values along ABD which belong to the 
lane wer region by 
\ particularly convenient formula for the circulation can now be given: 
rif P be any point on the stream-line BD and C is any curve which starts 
In Fig. 1 tl les of the rectangle EFGH are supposed to be remote from the aerofoil. 
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from P in the lower region, encircles the aerofoil and finishes at P in the 
upper region (without, of course, anywhere crossing the line A BD), then 


iT un > P ds = dp—dp, (9 


Cos 


cy 


é 
where the suffix P indicates values of ¢ and ¢’ at the point P. 
4. Since BD is a stream-line, the condition (8) must be satisfied at al] 
points along it and on both sides of it. Although the correct position of BD 
is unknown to begin with, and must emerge as part of the numerical solu- 
tion, nevertheless with any guessed position of BD, values of ¢ can be found 
by the usual processes of relaxation to satisfy equation (6) within the region 
ABDEFGHDBA together with the boundary condition (7) on the sides 
of the rectangle HF GH and the condition (8) on the double boundary A BD) 
There remains to be stated a criterion which provides a test, once such a 
¢-distribution has been found, of the accuracy of the curve chosen for BD 
and also a means whereby the initially guessed curve may be systematically 
modified, if necessary, so that this criterion eventually becomes satisfied 
When this has been achieved the solution to the flow problem is determined. 
The criterion which tests the accuracy of any guessed position of the 
stream-line BD is provided by the condition that the pressure along it 
should be the same when calculated from above (in terms of 4) as from 
below (in terms of ¢’). This condition implies that the velocities at any point 
of BD derived from ¢d-values and from ¢’-values should be the same; these 
velocities are both directed tangentially and, if s measures the arc-length 

along BD, the condition is that 
0h = 


Cs Os 


(10 


When this condition is satisfied it automatically implies, as should be the 
case, that (9) yields the same value for the circulation wherever P is taken 
to be on BD. 

The computational method of solution of this two-dimensional problem 
need only be summarized here, since it is similar to methods which have 
already been described (2, 3) for solving other problems in which there is 
an initially unknown boundary to the region of integration. We begin by 
making some guess for the position of the stream-line BD and then deter- 
mine the corresponding potential-distribution satisfying the governing 
equation (6) and the appropriate boundary conditions; from the values 0! 
the potential so obtained the quantities é4/és and é¢'/és are evaluated 
along BD, and then a modified curve BD is drawn on a basis of these values: 
thus, at points on the original curve BD where @¢/és is found to be smaller 
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than é¢’ és, it is appropriate (in Fig. 1) to draw the modified stream-line BD 


hove the original guessed position and vice versa. With this modified 
position for the curve BD, a correspondingly modified potential-distribution 
sdetermined and the computational cycle repeated until eventually a curve 
BD is found along which é¢/és and é¢’ /és are everywhere equal to each other. 
When this has been done the required solution to the flow problem is found; 
this is not difficult to achieve because the experience of a few trials quickly 
enables the right position of the stream-line BD to be determined. 
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5. Although other details of the relaxation method of solution of 
Laplace’s equation need not be described here,t the method of satisfying 
the boundary condition (8) along the line A BD deserves special mention 
because of certain features which are somewhat peculiar to the present 
problem. Fig. 2 represents a part of the line A BD together with a part of 
square relaxation net of side h: typical nodes of the net above AD are 
denoted by 0, 1, 2, and 3 as illustrated, and typical nodes below A D similarly 
v 0’, 1’, 3’, and 4’. P is the point on AD where it is crossed by the mesh- 
line 00’. The potential function ¢ exists only above AD where its nodal 


ilues are denoted by dp, ¢,, etc.; below AD the potential function is ¢’ and 


” 
tisdenoted at the nodes 0’, 1’, etc., by 45, ¢;, etc. At the node 0’, below AD, 
Is fictitious and its fictitious value there is denoted by ¢,,; similarly at the 


node 0, .bove AD, d’ is fictitious. 


They have been frequently explained elsewhere—e.g. by Southwell (4). 
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Laplace's equation it is 


F, = $,-+¢e+¢e+4y—4¢o> (1) 


the boundary condition that, at P, 


C d 


CV 


Q. (12 


If the slope of the curve AD at P is at an angle « to the 2-direction, down 
wards to the right, then (12) may be replaced by 
Ch od 


sin «+ — cos a 0, 
Cx cy 


for which we take the finite-difference approximation 


tan 
(4, $3) - a (do dy) 0. 


and thus the fictitious dy can be eliminated from (11) to yield 
tan a 


» 


F, d,-+4 dy ! ds 3dh5 | (4, ds) 


In a similar manner a residual formula may be found, in terms of ¢’-values 
at the node 0’: 

™ , , ae tan « 
Po = $1+4¢3+¢4—369—($1—93)- (14 


» , 


the numerical work can best be carried out, not directly in terms of the 


potential, ®, which is related with ¢ by the equation 


¢ (or d’) Ux+@® (or 0’); (15 


the presence of the aerofoil in effect induces a modification © which, wher 
added to Ux, yields the actual potential ¢ of the required flow. The formula 
(13) and (14) become, in terms of this induced potential, 


. ' tan a 
Fy = 0,4 0,+,— 30,4 Uhtana+(®,—®,) | 


1" , ’ ee : ,, tan a 
and Fy = 0+ 034 0,—30,—Uh tan a—(®,— 3) — ; 


where h is the mesh-length of the square relaxation net. In the numeric 


; mas tan a , ,, tana 
illustrations of this paper the terms (®,— 3) ——— and (®— 3) —— wer 


Now consider the residual formula at the node 0; in the usual way ¢y 


and contains the fictitious value ¢, which has to be eliminated by use of 


6. It became apparent, in the computation of the first example, that 


velocity-potential ¢ but instead, in terms of what has been called an induced 


Ux is the potential for the flow of an undisturbed stream of fluid so that 
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found to be so small that they had no effect on the last significant figure 
ioted; as a consequence the work could be carried out, and results quoted, 
n only one half of the complete field of flow, there being effectively anti- 
symmetry in the values of the potent ial about a line in the 2-direction drawn 
through the mid-point of the aerofoil. This circumstance is not, of course, 
feature of all problems of this kind (it would not, for instance, occur if the 
ilues of a happened to be larger than they are in the examples of this 
yper), and w hen it is not the case the solution has to be worked over the 
vhole field: this involves more numerical work to determine the solution 
but no greater difficulty of principles. 
7. Only one example has been worked in two dimensions, with results 
vhich are summarized in Fig. 3. The aerofoil was taken to be a thin plane 
finite rectangle represented in Fig. 1 by the straight line AB. Its angle 
f incidence to the direction of undisturbed flow is 2° and the values of the 
luced potential ® in Fig. 3 are given as multiples of 6000/U L, where L is 
lenocth of the chord AB of the aerofoil; these values of ® are anti- 
symmetrical (ef. section 6). Fig. 3 also shows the final position of the stream- 
ne BD, with values of « recorded at points along it. What was described 
section 4 as a process of drawing successively modified curves BD is, 
mputationally, a process of gradually changing these a-values. Fig. 3 


| on the finest net used; it is graded to a 


shows the solution obtaine 
smallest mesh-length of size h L/12 in the neighbourhood of the line 
{BD 
The solution was in fact obtained on three nets for which successively the 
nallest mesh-length was L/3, L/6, and L/12; a comparison of the values 
found on these successive nets enables an estimate to be made of the 
nvergence of the solution as the net-size is gradually reduced. Thus the 
lues of the circulation A, defined by (9), on these three nets were found 
be 0-058, 0-082, and 0-095 respectively, as multiples of UL; and the 
estimated value to which these successive results are converging as A is 


de smaller and smaller is 0-112U ZL. (1f the values found for AK on three 


successive nets are A,, A,, and Ky, then the estimated correct value for K 
sevaluated according to the assumption that the same ratio of errors on 
cessive nets is perpetuated as the nets are made finer and finer, i.e. that 
K,—-K K,—K' 
K,—K K,—K- 
\nother important quantity which can be evaluated directly from the 
tential distribution shown in Fig. 3 is the position, along the chord of 


erofoil, of the aerodynamic centre (or centre of pressure). If P is the 


cess pressure, above the uniform pressure far upstream, at a point in the 
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fow, and V is the velocity at this point, then Bernoulli’s theorem shows that 


where p is the fluid density. Hence at any point on the aerofoil the difference 
between the pressures, P and P’, on the upper and lower surfaces, is given 


hy 


p’_ p —Pp2_ yp’), (17) 


where V and V’ are the corresponding velocities. These velocities may be 
evaluated from the ¢-values given in Fig. 3, and so the pressure-difference 
17) can be determined at points along the aerofoil section AB (Fig. 1):+ 
the position of the aerodynamic centre has been found by numerical integra- 
tion (using Simpson’s rule) along the chord AB; it is at a point distant 
):256.L from the leading edge. 

8. The example in the previous section can be solved exactly, so that we 
ire able to test the accuracy of the relaxational solution. Thus the correct 

ilue for the circulation K is 0-110U L, and the value obtained by relaxation 
methods is therefore seen to be in error by about 2 per cent. Similarly the 
orrect position of the aerodynamic centre is at a point distant 0-250L from 
2 the leading edge; the relaxational value is therefore in error by 0-6 per cent. 
5 of the length of the chord. 

The accuracy obtained by relaxation methods is considered to be very 
satisfactory, and this example suggests that similar results obtained by 
the method can be accepted for other aerofoils where the shape does not 
permit of an exact solution (for although the aerofoil in the example of the 
previous section is of special shape from the orthodox analytical point 
if view, no shape has any particular significance from the point of view of 
relaxation). We give here no further examples of two-dimensional flow, but 
turn instead to consider problems in which the flow is essentially three- 
limensional. In such cases (when the aerofoil is of finite span) exact solutions 
in hardly be contemplated, and any theoretical assessment of aerodynamic 
properties must rely on some form of approximate or numerical method. 
The satisfactory nature of the two-dimensional solution obtained in section 7 
by relaxation methods, together with the fact that these methods can now 
be applied in three dimensions, would seem to imply that they provide a 
reliable means of dealing with the flow round three-dimensional aerofoils. 
This expectation will be seen later (in section 13) to be confirmed. 

9. The extension of the method into three dimensions can be described 
with reference to Fig. 1. The rectangular boundary EF FGH is replaced in 
three dimensions by a rectangular parallelepiped, all of whose faces are 


At A and B (17) is not used, the pressure-difference being taken to be zero. 
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remote from the aerofoil; the line AB in Fig. 1, representing the aerofojj 
becomes in three dimensions a fixed (and finite) internal surface, and thy 
stream-line BD becomes a surface (or sheet) extending downstream fro; 
the trailing edge of the aerofoil. The position of this sheet is not know 
initially but has instead to be found as part of the computational solutio; 
in a similar way to that in which the line BD has to be found in ty 
dimensions. 

In three dimensions we choose orthogonal axes of x, y, and z with the x-axis 
parallel to the direction of undisturbed flow, the y-axis perpendicular to th 
span of the aerofoil, and the z-axis parallel to the span; the equation gover 
ing the velocity potential ¢ is then 
2d, Od 4 Og _&¢ 7 Y (18 

C 


7 a? 
x oz" 


9 | ' 


oy” 
and the boundary conditions to be satisfied on the faces of the parallelepipe: 
(the external boundary) are that, on the two faces which are perpendicular 
to the direction of undisturbed flow, 


Ch . 
ft UE (19 
Cx 
and, on the remaining four faces, 
od od 
— 0 or -= Q. (20 
Cy Oz 


The condition to be satisfied on the surface of the aerofoil (the interna 
boundary) and on both sides of the sheet is 


0. (2] 


where @/év denotes differentiation along the normal at any point of the 
aerofoil surface or sheet. For any chosen position of the sheet a potential 
distribution, ¢, can be found satisfying equation (18) and the conditions 
(19), (20), and (21). Before stating the second condition which has to be 
satisfied on the sheet in addition to (21), and which is used in order that the 
correct position of the sheet may be eventually determined, we consider first 
how (21) may be applied for the determination of 4 with any guessed 
position of the sheet. 

10. Fig. 4 represents the three-dimensional counterpart of Fig. 2: pat! 
of a relaxation lattice is shown in the neighbourhood of the sheet whichis 
represented by the two lines APD and JPK. Points of the lattice are 
denoted by 0, 1, 2, 3, I, and II (above the sheet) where corresponding 
potential-values are dy, ¢,, etc., and by 0’, 1’, 3’, 4’, I’, and II’ (below the 
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n the sheet where it is crossed by the lattice-line 00’. 


2 




















Ri 
— = ie) I 
(A | 
—_—_—> X 
( \P (J) 
3) | (D) 
z<" - ie si sj m, 
3 —T 1 
| 
|, 
Fic. 4. 
In the usual way for Laplace’s equation (18), the residual formula at 
oOlIntT V 
F, Pi 7 P2 ps do py T by bho, (22) 


ntaining the fictitious value which has to be eliminated by means of 


the condition (21) holding at P. If the (downward) slope of the curve AD 


tothe x-direction is at an angle « and the (downward) slope of the curve JK 

the z-direction is at an angle £, it follows that the direction ratios of the 
normal to the sheet at P are (tana, 1, tanf) so that at P the condition 
2] may be re placed by 


CO c ad Cc dh ) 
tana tanp = V, 


C2 cy Cz 
lor which we take the finite-difference approximation 


tan B 


» » 


tan a 
d f 
} ) Py 


(dy py) 0, 


MD 
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and thus the fictitious 4, can be eliminated from (22) to yield 


tan a tan B 


+ (dy dy4)- ——- ie 


Ky = $)+$2+¢3+¢,+b11—5h9+ (4,—9s) 





As in two dimensions the numerical work is best carried out in terms of ay 
induced potential ® which is again related with ¢ by equation (15); in terms 
of ® the formula (23) becomes 


F, = 0,4 0,+ 0,4 0,+ 0,,—5®,+ Uhtan a4 


tan a te in f 


-(®, %,)—=-+ (®,—®,,) (24 


and in a similar manner a residual formula may be found, in terms of 
®’-values, at the node 0’: 


Vv , , , ’ , ' , — , y 
Fy, = 9,4 0,4 0,4 0,4 0,,—5®,— Uh tan a— 





tan « =— 


(D; ~',) - — — (9; —; ) 
For the reason already given in section 6 there is effectively antisymmetr 
in the values of the potential—here about the za-plane drawn through the 
mid-point of the aerofoil; as the aerofoil is itself symmetrical, there is also 
symmetry in the values of the potential about the zy-plane drawn through 
the centre of the aerofoil. As a consequence the computations in ow 
examples could be effected, and results quoted, over only one-quarter of 
the complete volume of flow. 

11. To ensure that eventually the correct position of the sheet is found 
we must satisfy, in addition to (21), the condition that the same value of the 
velocity at any point on it must be obtained from the ¢-values above as 
from the ¢’-values below (the velocity, in the precise meaning of the word 
need not be the same above as below since its direction is not necessaril) 
continuous across the sheet, although, of course, both above and below itis 
tangential to the sheet). The velocity, V, at a point P (Fig. 4), above the 


sheet, is given by = (24 2 (9 2 (28 P 
Ox} ° sé ' Naz]? 


with a similar expression for the velocity V’, below the sheet, in terms of ¢. 
The values of V and V’ can be found at points on the sheet, from any 
relaxational solution for the potential obtained using any guessed position 
of the sheet; this position can then be modified gradually to correct an) 
discrepancies between pairs of values of V and V’ in a similar manner t0 
that described in section 4 for the gradual modification of the position of 
the stream-line in two dimensions. 
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12. We specify the position of the sheet numerically by assigning values 
othe distances typified by 0’ P (Fig. 4): thus 7) measures the distance 0’ P 
;a fraction of a lattice-length, and 7,, 73, 1, and yy, measure similar 
iistances as Shown. Then, when required, as in (24) and (25), we may use 


the finite-difference approximations: 


3 1 Ki % mu 


») 


tan a and tan - 








Fie. 5. 


13. Two examples have been worked in three dimensions: both are for 


flow round thin (plane) aerofoils. The first case deals with an aerofoil of 


rectangular profile, with aspect ratio four, placed at an angle of incidence 


to the direction of undisturbed flow (this aerofoil is thus similarly 


tuated to the two-dimensional aerofoil discussed in section 7, and differs 


only in having a finite instead of an infinite span). The second case deals 


with a swept-back aerofoil, of aspect ratio three, whose plan is shown in 


Fig. 5; the angle of sweep is 45°, and the angle of incidence again chosen 


» 


Details of the relaxational solutions, and results in the form of potential 


ues at the points of a relaxation lattice, are given in section 15. We 


irst quote some important results derived and compare them with those 


uned by Falkner (5, 6) by use of an entirely different (‘vortex-lattice’ ) 
ipproximate method. 
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The circulation round a closed curve C is again given by the contoy 
integral (9), and interest is attached to the variation, along the span, of th; 
circulation round sections of the aerofoil by planes perpendicular to the 
spanwise direction. 

In both cases the relaxational solution was obtained at the points of thre 
lattices for which the smallest mesh-length was L/3, 1/6, and L/12 succes. 
sively, where L again is the length of the chord of the aerofoil (cf. section 7 
The estimated values, to which the circulation is converging at sections 
along the span, are shown in Table I; the values of K are given as multiple 
of UL, and position along the span, from its centre, is measured non 
dimensionally by z’ (= z/S, where S is the length of the half-span). 


TABLE I 


Rectangular aerofoil Swept-back aerofoil 
ispect ratio four aspect ratio three 
2° Kk ic. 2° K ALC. 
° 0*070 0°245 ° 0°045 0°34I 
1/24 0"070 0°245 1/18 0°040 0°305 
2/24 0°070 0°245 2/15 0°047 0°255 
3/24 0070 0°244 3/18 0048 0°273 
4/24 0070 0°244 4/18 00458 0° 264 
5/24 0°069 0245 5/15 0°049 0°257 
6/2 0009 o°244 6/18 0°049 0°245 
7/24 0°065 0°242 7/15 0050 0°245 
8/24 0°007 0°243 8/15 o°050 0°239 
9/24 0066 o°24I g/1d 0°049 0°230 
10/24 07005 0°239 10/15 0°049 0°232 
11/24 0064 0°242 11/18 0°045 0°222 
12/24 07003 o*24I 12/15 0°040 0*205 
13/24 0061 0°237 13/15 0044 0195 
14/24 0060 0°234 14/15 0040 Or175 
15/24 07055 O23! 15/15 0°035 or155 
16/24 0055 0°227 16/15 0°029 0°134 
17/2 0°053 0225 17/15 0°020 o°105 
15/24 0050 0°223 15/15 0*000 
19/24 0°047 o°221 
20/24 0°042 0°217 
21/24 0°037 oO°213 
22/2 0°030 0°209 
23/24 0°021 0°203 


24/24 0*000 





For comparison with Falkner’s values, these results are presented graphi- 
cally in Figs. 6 and 7; in both cases the normalized circulation, A’, is plotted 
against z’. The normalized circulation is given in terms of K by the formula 

1 


Te 
Kk’ x | K dz’. 
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Kia. 9. 


Figs. 8 and 9 show how the position of the local aerodynamic centre varies 
along the span; in both cases this position is measured by its fractional 
distance} along the chord from the leading edge, and a comparison may be 
seen with Falkner’s results. The close agreement between the two methods 
to be seen in all of Figs. 6-9 lends confirmation to our belief that the 


+ Values of these fractional distances are also recorded in Table I in the columns headed 


‘A.C.’ The quantities recorded in Table I were evaluated for each section of the aerofoil 
in the same manner as that described in section 7 for the two-dimensional case. 
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relaxation method yields results in three dimensions which may be accepted 
to be of the same order of accuracy as those which it is known to yield in 
two dimensions 
14. Before we discuss the actual relaxational solution for ¢ at the points 
fa cubical lattice, an account should here be given of how such a lattice 
nay be graded from a region where its side is of length h to a region where 
ts side is h/2; this device of graded lattices in three dimensions is similar 
to the device of graded nets in two dimensions which has been described 
\llen and Dennis (7). Fig. 10 illustrates a small portion of the lattice 


the region where the coarse part is connected up with the finer part: 





Fic. 10. 


the coarse lattice extends from the right in Fig. 10 up to the yz-plane FEK L, 
vhilst the finer lattice extends from the left up to the yz-plane CDJH. 

order that the lattice may be graded from the coarse to the finer size 
we must be able to define residuals at all lattice-points which are introduced, 
nd, in Fig. 10, residuals can be immediately calculated at all points in the 
plane FEK L according to the usual formula 


Fy = Sw 


a > “1 — 6u, (26) 
where § w, stands for the sum of six w-values at points distant h from the 
typical point 0. Similarly, residuals may be calculated according to the 
same formula at all points in the plane CDJH: at the points in that plane 
typified by C, D, H, and J, ¥ w, again stands for the sum of six w-values 


it points at distance h; but at the points typified by Q, R, S, 7, and U, 
>. 


; Must stand for the sum of six w-values at distance h/2. This requires 


5092.21 


H 
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the introduction of lattice-points situated like B, G, M, N, and P, inte 
mediate in position between the coarse and fine parts of the lattice. |; 
remains therefore to state formulae whereby residuals may be calculate; 
at these intermediate points; they are of two kinds, those like G which ar 
placed at the centres of cubes of the coarse lattice, and those like B, My 
and P which are at the centres of faces of such cubes. We shall demonstrat, 
how residuals may be calculated at these points by considering one of eac} 
kind, i.e. the points G and B. 

So far as B is concerned we have at once the usual finite-differcog 


approximations: 


Cw wW,+W;—2Wz 
cy” (h/2)? 
and, in the plane CDEF, 
. ‘ 
Ow  0o-w Wy Wytwy, Wp 4, 
cz? a? (h/v2)? : 


Adding these two results together, we find that Laplace’s equation 


satisfied if P 


=(W 4 Wy) + UY WytUp_t uy SWp 0: 

hence the residual at B is defined by the formula: 
’ » 9 ; 9 ! 9 9 5 Pass i 
F, 2(Wy W,,) WotrwWp Wy Wy, SW. Zz 


At the point G, it can easily be shown, by means of an expansion i1 
triple Maclaurin series, that 


iw iw u“ 


( D u 


7 Wry Wy FW Ew, Sw, h?(V2w),, | O(h4): 


thus, with neglect of the terms O(h*), we have the finite-difference approx 


mation to (V*w),, leading to the definition of a residual at G by the formul 


KF uw 7 uw 


. . ‘ . . 2 . a 
0 ( prUp rept Uy t+ Wy UK + Wy — 8G . 


as 

Hence, by means of formulae of the types (26), (27), and (28) residuals 
can be calculated at all points of a graded lattice. 

15. A number of such gradings of the lattice were made in the solutions 
for the velocity-potential. The complete solutions are too extensive to bi 
given here (for the rectangular aerofoil the finest lattice contained ov 
1,900 points, and for the swept-back aerofoil over 1,750 points). Instea’ 
we only show, by way of illustration, in Fig. 11, a part of the finest lattice 
in the immediate neighbourhood of the swept-back aerofoil. The values 
recorded are those of ®, and a single grading of the lattice can be seen rout 
the edge of the diagram. 

Only one point of detail in respect of the solutions calls for special mentio! 


here and this concerns the spanwise extent of the sheet downstream of tht 
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erofoil. We have, in section 11, described how the sheet may be adjusted 

its correct position, but we were only dealing there with the gradual 
letermination of the position of points of the sheet along lattice-lines in the 
direction (as the sheet is throughout inclined at only a small angle to 

pnlanes its position is primarily determined by knowledge of the position 
f such points). In plan view the sheet extends downstream from the 
railing edge of the aerofoil in the positive x-direction indefinitely (as does 


| 


stream-line in two dimensions); but, also in plan view, the extent of 
the sheet in the spanwise direction is not necessarily constant but instead 
iy be expected to spread slightly, downstream from the aerofoil. This 
need not cause any extra difficulty because any assumption of the 
spanwise extent of the sheet, greater than the actual extent, does not 
produce any error in the solution; the assumption that the sheet extends 
eyond its actual limits only allows for a discontinuity where in fact none 
xists (it does not compel a discontinuity to be present). In the two solutions 
vhich have been given here for flow round three-dimensional aerofoils the 
sheets were not found to extend in the spanwise direction as much as one 
ttice-leneth 
We have not made any allowance for instability of the sheet, judging that 
tendency to ‘roll-up’ will only become appreciable at a distance which 
s sufficiently far downstream to justify its neglect. 
16. In conclusion we would emphasize that the method described in this 
per is the one which it would be necessary to adopt for the determination 
f the flow round a three-dimensional aerofoil in the most general case. 
ertain modifications could have been made, in the particular examples 


r which results have been quoted, which would almost certainly have led 


to simplification of the computations involved. Thus it would be better 
hen the aerofoil is thin and plane to orientate the relaxation lattice so that 
lattice-lines in the 2-direction are chosen to be parallel to the chord of 
he aerofoil instead of being parallel to the direction of the undisturbed 
w: in general when the aerofoil is thick (and not symmetrical in section) 
ttice is better orientated with the x-direction chosen to be parallel 
the direction of the undisturbed flow (as shown in our examples). 
\nother, and more considerable, shortening of the work involved, could 
made if one is prepared to accept as a justifiable assumption ab initio, 
nany o1 ane perpendicular to the span the circulation has the same 
iue round all curves enclosing the aerofoil; we did not make this assump 
n, but our results do verify its validity for the two examples which have 
en worked 
Grateful acknowledgement is made to the D.S.I.R. for a grant which 


enabled this investigation to be completed. The paper describes in part the 
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THE DIFFRACTION OF A DIPOLE FIELD BY A 
PERFECTLY CONDUCTING HALF-PLANE 


By T. B. A. SENIOR (Cavendish Laboratory, Cambridge) 





Ree d 1 November 1951] 


UMMARY 


g a method due to Bromwich the diffraction of a dipole field by a perfectly 
lucting half-plane may be expressed in terms of the known solutions of a scalar 
em, but only for a dipole whose axis lies parallel to the diffracting edge. An 

ternative method is developed which is valid for all positions of the dipole and 

m the resolution of the dipole field into plane waves whose diffraction by the 

may be independently studied. The particular case considered is that of an 
tric dipole with its axis normal to the half-plane. 


1. Introduction 
[ue diffraction of a dipole field by an infinitely thin, perfectly conducting 
ilf-plane is a problem which, though fundamental in diffraction theory, 
isnot been adequately treated. There is, however, a corresponding scalar 
roblem, involving solutions of the wave equation which are finite and 
mtinuous everywhere except at a given point, where they behave as 
*'’R for small R, and which vanish, or have vanishing normal deriva- 
tives, on the sheet. Solutions have been given by many authors, notably 
slaw (1), (2) and Macdonald (3). 
The half-plane problem is clearly a special case of diffraction by a wedge. 
fo this also there corresponds a scalar problem, solutions of which have 
een used by Bromwich (4) to deduce the field of a dipole in the presence 
1 wedge. Essentially his method consists in choosing a Hertz vector 

I= (/U,mU,nU) where 1, m, n are the direction cosines of the dipole and 
is a solution of the scalar problem, but only if the dipole has its axis 
rallel to the edge can the boundary conditions be satisfied. 

It seems feasible to apply this procedure to diffraction by a half-plane 
nd no restrictions upon the orientation of the dipole are then necessary. 
{san example, consider the particular case of an electric dipole whose axis 

normal to the sheet and parallel to the y-axis of some coordinate system. 
hoosing the Hertz vector II (0, U,0), the boundary conditions on the 
sheet are satisfied if el’/ey vanishes there. The diffracted form of U, 
owever, has a singularity of order r! at the edge of the sheet, giving rise 
singularities exceeding r—! in the field components. Recent work on edge 
singularities by ( opson (5), Meixner (6), and Jones (7) indicates that such 

solution of the vector problem is inadmissible. Moreover, use of this 
Hertz vector presupposes the absence of any component of the magnetic 


(Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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field in the y-direction. Whilst this is true of the incident field, it is no 
necessarily true of the field scattered by the sheet, and a component H, wil 
be shown to exist. 

Nevertheless, the method can be applied to diffraction by a half-plan 
if the dipole lies parallel to the edge, the solution so obtained satisfying al 
the boundary conditions including those at the edge. This may be attri. 
buted to the fact that the incident and scattered fields then have the sany 
polarization. Thus, for an electric dipole parallel to the edge, the incident 
magnetic field has no component in that direction and a scattered field of 
the same form, generated by currents in the sheet flowing parallel to th 
edge, is consequently a possible solution. For a magnetic dipole a simila 
result follows by invoking Babinet’s principle. But for other positions of 
the dipole the induced currents are such as to generate in the scattered field 
a component which, in the incident field, is zero; the extension of Brom 
wich’s method then breaks down. 

An alternative method is therefore developed which relies upon th 
expression of the dipole field as a double integral over plane waves whos 
diffraction by the sheet can be determined. The case considered is that of 
an electric dipole with axis normal to the sheet, but the method is equall) 
applicable to a dipole orientated in any way. The analysis given refers ti 
the component of E normal to the sheet; the other components may b 
obtained in a similar manner, only their final expressions being presented 
here. 


2. The method of analysis 

Consider an infinitely thin, perfectly conducting sheet occupying the half 
plane y 0,2 0 of some cartesian coordinate system (x, y,z). An electri 
dipole whose axis is normal to the sheet is situated at the point (Xp, 4,2 
in terms of the Hertz vector II (0, 11,,,0) where IT, (e*R ER), R being 


the distance from the dipole, the free-space field of the dipole is given by 


. eu, ei eq ae ell 
E v, —#+Il,, 34 Hm = ski yr, 0, ——*3,.. @ 
cxcy cy? ” Gydz oz Cx 
where Y 1/Z is the intrinsic admittance of free space. The m.k.s. systen 


of units is employed and a time factor e’* suppressed throughout. 

We first seek a representation of the field (1) in terms of plane waves 
As is well known, 
ikR ae . 


é » : . ) 
7 | do | eikR{cos(a—O)cos @ cos B-sin Dsin BJ egg BdB, (2 
kR 2m 


Dp ¢ 
. . 


7 ix 


where R, ©, and ® are spherical polar coordinates with the dipole as origi" 
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For our purpose, however, we require a modification of (2) in which the 
ntegrations are over ‘steepest descent’ paths. This we obtain as follows. 


wince 


etkR cos(a—O)cos ® cos B dy J,(k R cos ® cos B) 


mation (2 becomes 


J,(kR cos ® cos B)e-** 51 ® Sin B eos B dB 


H\’(k R cos ® cos 8) + H\?(—k R cos D cos B) 


e—tkR sin ® sin B agg B dB 


H' (kR cos ® cos B)e-t*# sin @ sin B eos B dB, 
») 


I 
where P(47) denotes the steepest descent path passing through the angle 
7. Inserting Sommerfeld’s integral representation of the Hankel function, 
ve find 
| | , Ricos(a—@)cos ® cos B—sin ® sin B cos B d vd 8. 
A R “71 . « 
I 5=f 


nd the integrals are now of the required form. In terms of the original 


ordinates eile mite 


( sacos P+y sina cos P—zsin b) F'(y B)eos BdadB, 
Pum) B=Pus (3) 
vhere FY x. B) , } COS a COS 6+Y_ SiN a cos b—z, sin 6 


\ssuming that differentiation with respect to a, y, and z under the signs of 
ntegration is permissible, we are thus led to the following expressions for 
the components of the dipole field 


. 


E | sin x COs a cos?8, 1 —sin?a cos?8. sin asin 8 cos B} 
Niky) 
i P T 
é x cos P+y Sin a cos b—z8IN B) F'(y B)cos B dadB, 
, ; (4) 
H = | Y sin 8, 0, —Y cos «cos f! 
a 


) 


- RR ; R , o 
etka cos P+y sin a cos P—26in B) H'(q~ B)cos B dadB. 
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Equation (4) represents the dipole field as a double integral over three. 
dimensional plane waves whose diffraction by the sheet can be determined. 
It is convenient to reduce the diffraction problem to one involving two. 
dimensional, plane-polarized fields. 

Now any two-dimensional solution of the wave equation gives rise toa 
three-dimensional solution on replacing k by kcos8 and multiplying by 
e-tkesinB Tf V is such a three-dimensional solution, an elect romagnetic field 
may be derived from it by taking V as the z component of an electric o 
magnetic Hertz vector whose x and y components are zero. In this way we 
obtain two fundamental types of field: 

(i) E-polarized in which H, = 0, 

EO (Spee isin B eV 


: Veos*B}, H®) 24 : vo ov 0): 
k Oe kk wy a 


k Cy ex 
(ii) H-polarized in which E. = 0, 
isin B eV isin B eV 


k ox kay 


— a eV eV 


, 0}, H® 
| 


! . Vcos*g). 
CY Cx 


(6 
Consider that solution of the wave equation which corresponds, in the 
| | 
manner described above, to eos x+y sina) Tt ig eik(rcosacos B+ysinacos B—zsinp, 
in terms of which 
E® (—cos asin 8 cos 8, —sin «sin B cos 8, —cos?B) » 
v pik( cos a cos B+y sin a cos B—z sin B) 
H® = (Y sinacos B, —Y cos «cos, 0)e**@ cos « cos B+y sin a cos B—z sin B) 
E®) (Z sin «cos 8, — Z cos «cos B, O)e*@ 08 «C08 B+y sin a cos B—z sin B). 
H® = (cos asin cos B, sin asin B cos B, cos?8) » (8) 


“a ik(a cos « cos B+y sin a cos B—z sin B) 


When 8 = 0 these reduce respectively to E- and H-polarized fields in two 
dimensions. This suggests that the solutions of the diffraction problem fo1 
the incident fields (7) and (8) may be deduced from those in the case B = 0), 
a procedure which is certainly valid if the diffracting structure is two- 
dimensional. Thus, in the two-dimensional solution for Z_, replace k by 
k cos 8 and multiply by e-‘**si"®, When substituted for V in (5), this forms 
the solution of the diffraction problem for the incident field (7). And 
similarly for H-polarization. 

In order to combine the fields (7) and (8) so as to produce one whose 
magnetic vector has no component normal to the sheet, constants A and 
must be chosen such that 


H, = \H0+pH® = 0. 
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This condition is satisfied by A = sinasinf and pz Y cosa, and the 
resulting field is then 
E = {sin a cos « cos*8, —cos B(1 —sin?a cos?8), —sin «sin B cos?B} x 


v pik(x cos « cos B+y sin a cos B—z sin B) (9) 


. > : 
4 > 4 » > COS ‘os Bays s B-—z 
H } sin BcosB. 0. } cos v COS*B eth co x cos B+y sin a cos B z sin B) 


ik? 


ONES 
“TI 


If these expressions are multiplied by | )F(a.B) and integrated with 


respect to « and 8 over the path P(}z), we recover the field of the vertical 


electric dipole. Lt follows that 


E, H) - : | sin asin B(E™, H®)-+ 
- Pim) B=P(in 
Y cos a(E®, H™)} F(a, 8) dadB (10) 
lin particular 
*L.3 ” ° 
E, = | sin asin BE) +-Y cos vE\?)} F(a, B) dadB. (11) 
= . " > 


This forms the required resolution into plane waves. 


3. The diffracted field of the dipole 

Equation (11) serves to reduce the solution of our problem to quadratures 
rovided we can determine the diffraction effect of the sheet upon the 
elds (7) and (8). For this purpose let us return to the study of two- 
limensional plane polarized fields. 


Consider an #-polarized wave incident upon the sheet and let 
E pik(a cos a+y sin a). 
the affix “7° will be used to denote the incident field. The solution of this 


liffraction problem is well known, the most convenient form here being 


that given by Clemmow (8) 


. a 1 .. - 
E, = Bi—~ | — 2° 2 3Y ¢-aiecosysivisin y) dy (12) 
qT COS «& COS Vv 
/ 


( 


where the path C is as shown in Fig. 1. Generalization of (12) as indicated 
»4 


nsection 2 then leads immediately to the solution for the incident field (7). 


Similarly, from a two-dimensional H-polarized wave in which 


H ( k(a2 COS a+y Sin a) 
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and whose diffraction is given by 


© oes ne la 
Cos 5) whl i # ik(a cos y+|y|sin y) 


H, = Hi4* 9 


a \y| | cosa+cosy 
F 


dy (13 


/ 


(Clemmow, 8), we may deduce the solution for the incident field (8). 











C 
A 
0 = J, 
> ‘ 
a Pia) } 
c 
Fic. 1. 


Inserting these results into equation (11) we arrive at the following 
expression for the y-component of the electric field of the dipole in th 


presence of the conducting sheet: 


2 . 2 
E, i | | (1 —sin?« cos?8) >» 
= x= P(47) B PGt) 
v eik(e cos «cos B+y sin x cos B—z sin B) F( xv. B)cos B dadB 
je2 y r » re 
2x2 y| | | | 
x= P(r) B=PUm) y=" 


PERO Pe Rr TO ee, tae a 
{COS 5a COS sy COS a COS y—SIN 5aSIN 5y SIN asin y sin*p} 
| COS a-4 COS Y | 

x e—ik(e cos y cos B+\y\sin y cos B+2 sin B) F( x,B)cosB dadBdy. (l4 
We must now consider the evaluation of these integrals. In the ensuing 
analysis we shall require the polar coordinates (7, @) defined by x = rcos¢ 
y = rsin@; as before, a suffix ‘0’ will signify the coordinates of the dipole 
Assuming that differentiation under the integral signs is permissible, tht 
double integral in (14) reduces to 


/ 22 \ e-ikR 
( a | ie (15 
CYCY, kR 
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whilst the triple integral becomes 


l 


{sin } L, 


. 2 
YSIN sy o~ 


> | | | 
= ; ; (cos a+ cos y CYOYy 


Pan) B= Pam) y= 
ee, Pe ee em a? sin v\ k2 
ad tas 5y COS ~ COS Y sin 5a SIN 5ysmM asin y Y\ 
cos a+ cos y ly| J 


e—tk(x cos y COS f y sin y cos B+2sin B) F'(y, B)cos B dadBdy. (16) 


For separation into a geometrical optics-field and a diffraction field, C must 
be distorted into a steepest-descent path. The saddle point for the region 


0 is 3 6 and for y < 0, 4 2z—6. It may be shown that the pole 
7—a is captured (in a positive sense) when 7—a lies to the right of P(@) 
y > 0) or to the right of P(27—@) (y < 0). Since a lies on P(4z), the 


residue provides the following contribution to the geometrical optics-field : 


| . ke): i for 0 < 4n, 


(17) 
ikR 


ou , € . . 
k?) , for 0> 3x, 
CYUCYs kR iy 


vhere S is the distance from the image of the dipole in the plane y = 0. 


Consider further the triple integral in (16), where now the y-path of 
ntegration is P(@) or P(27—8@) according as y > 0 or y < 0 respectively. 
Regarded as a function of «, the integrand possesses a pole «a am—vy at 


vhich the I'é sidue IS 


| c- r ikS 
k2) —, forge (0, 
2mr1\cycy, kS 
(18) 
| e-tkR 
ke) , tory < 0. 
2m1\cy fo kR 
Deformation of the path of integration so as to pass through the a-saddle 
int (a #,) may require this pole to be taken into account; it is found 


it the pole must be included in a positive sense when 
= cs ) 5 one 3. 
XL 4 6, o1 a A, G 71 


nd ina ne iTl1VE sense when 


7—O iy lyr or 37 <6 a7+6,. 


0 pa 
[he residues in these cases provide additional contributions to the geo- 
letrical optics-terms which finally become 


¢ kh e2 e—tks 


| , k2) | f:2 —, for0<0< r—6h, 
CYCY, | kR \CYCY, ks 
4 : \¢ kR ( 19) 


\ouc Yo 
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The remaining integral represents the diffracted field. Denoting it by J 


we have, after some reduction, 


Sa? \cos 3 (a ryt 6-4 5) 


s= P(0) 


Dh : = | } | | Cos(a+-y- 6+-45) 


Y 
/ 


cos(a—y—86-+-6,) 
4 € 
cos $(a—y—6+-6,)} 


ik{r cos y cos Bry cos a cos B+(z—z9)sin B cos B d yd Bdy 


l a9 . > . 


zs fsec A(a+y+0+6,)- 
stay, | | | Seeaty+0+6) 


. © 
x= P(0) B= P(47) y= P(0) 


ikire y cos B41, COS “os B4(z—z2,)sin B} 
sec 3(a y A 6.) te ik{r cos y « Bir, Cos x cos B4(z—29  P} cos B dadBdy, 


This may be more conveniently written as 


kh 
i or {A(0,)+A(—8@,)} 
l 2 | a2 
: | = Kel B(99) — = / : 1 pel B(—4,) — (20 
2kileyey, j 2rkileycyy j 
where 
Ik i r . 
A(@,) ami | | COS ( oe me a 6 Ay) e 
7 x= P(0) B PGn) yY *P(O) 
Z e—ik{r cos y cos B+r9 cos a cos B+(z—z9)sin B} cosBd vd Bay 
and 
k ‘ . . 
B(@,) om | | | sec 3( X-+ Y t Q 4.) . 


. . 
x= P(0) B= Pin) y= P(0) 


v o—ik{r cos y cos B+rg cos x cos B+(z zosin B} egg B dadBdy 
f ply. 


4. The reduction of the solution 

We first reduce the integrals A(@,) and B(6,) to forms analogous to those 
given by Carslaw (1) in his discussion of the corresponding scalar problem. 
Consider the integral A(6,). Make the formal change of variable 


(a, B,y) > (a,B,5) where 6 = y+a. 


In order to determine the path traversed by 5, choose some fixed value of «. 


Since y lies on P(0), it follows that 6 will lie on a similar path which, though 
displaced by an amount «, can be shifted back into the original steepest 
descent path. Similarly for each value of « on P(0). Hence 
k ; re ; ; 

| | cos $(6+-0—8@,) » 


A(6,) 


x= P(0) B=P(ir) 6=P(0) 


z e—ikir cos(S—a)cos B +19 cos « cos B+(z—z9)sin B} cosBd yd Bdd. 
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The a- and 8-integrations can now be carried out immediately to give 


, o—tk 9) +2rre COS 0 : _ 
1(6,.) — . cos $5 cos $(8@—6,) dd. (21) 
J qr?+ro+ (z—z29)*+ 2rrg cos 0} . = 
P( 


Nin ilarly we obtain 


. ik 2rre cos 8 = = 
é 0 cos 40 cos 3(6—6,) ,. 
B(6,) -- 2 0 5 (22) 
rere 2— 2 9)*+-2rn cosdo} cos d+-cos(6—86,) 


this form A(@,) is identical with Carslaw’s basic solution of the scalar 
blem, apart from a constant factor. 
Further simplification of both A(@)) and B(@) is possible using analysis 
sed on that of Macdonald (3). Thus, consider A(6,). Since 
by . 
s 


» y 
é 


8 2—U~/ cs 
U (27) Ss? 


(23) 
here a is a real positive quantity, A(#,) may be written 


1(6 exp| s/2— (k?/2s){r?-+-r2 + (z—z,)*+ 2rr, cos 8}] > 
\(47) J, 


ds 
< cos $6 cos }(8—6,) dd a (24) 
si 


Now perform the integration with respect to 6. Using the substitution 


T V2 e-*7/4 gin ae) (25) 
C08 8 agg 15 ds V2 eitl4ek*/s)rro | eik®/syrrot® dr 
P(0 os 
’ ar 
‘ =) (k*/s)rro 
k AN rl 
Hence 
1(@ 1/A ) [e/9 k2 Do\S lp - \2 ~ ~ \2) ds 
A(G,) i cos $(86—6,) exp| 8/2 —(k*/2s)}(r-+19)°+(z 2o)"}]- — 
8\/(779) 
cos 1(0—6,)H(k R,), (26) 
'o) * 
where R? 1)? + (z—29)? 
In like manner B(@,) may be written 
) f c j 9 9 6 . i 
B(4,) - | exp| s/2—(k?/2s){r?+-r5+ (z—z»)?+- 277% cos 8}] > 
N\ST) J 
PO 0 


cos $6 cos }(8—6,) a®. (27) 


’ cos 8+ cos(0 Ay) si 
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The substitution (25) then gives 


P 2ng 1S ang 1 if] ' (k?/s)rr,r° 
2 >: COS 50 COS ¢ 0 . s , ellk*|s)rrg7 
e—(k*/s)rrg cos 0 2 2 0) ae e—it/4o (k* ST 0¢ | - dr 


. cos 6-+-cos(6—8@,) 7? te? 
P(0) 
| [Da . 
14 baa e(k?/s)rrg cos 45) | P A*/28 q) 
NV \s | 
€v( 2779) 
where e« v2 cos }(@—@,) and the upper or lower sign is chosen according 
as ba < arge < 4m or ba < arge < 37, respectively. Inserting into 
(27), f a ran 
° i i l 9 9 9 ds 
BiG) = tik | | explis——aerepr! a® 
| 28 | s* 
€ex(27rr9) 0 
C HOE /(k2R2+A?2)} 
& | 2 ay 
b \ (k2.R?+-A?) 
€4(27r79) 
Lark HP(kReoshp) du + for 6 < 7+, (28) 
LR 
where Lp = sinh 1/2070) cog 3(0 ,)|. 
, | R ” 
Similarly 
B(—4,) Lark | H?(kScoshp)du + for 6< 7—@, (29) 
‘Ls 
ae 
where Ls sinh Hew 0) cos 1(Q4 4,)}. 
< 2 


[t only remains to introduce these results into equation (20). Since 
Di ¢ iku 


| H'?(ku cosh p) du — 


the changes in sign of B(@,) and B(—@,) may be conveniently combined 
with the geometrical optics-terms (19) to give the following expression for 
the total field: 
y ik 2) » 
E, cos $4 cos $6, H\?(k R,)-4 
\ (779) 
al oO = b 19)/7, D 
—k H(kReosh p) du 
1 I I 
2\CycYy 7 
PR 


va 


( <l + ke) H(kS cosh p) dy. (30 
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It is interesting to compare this solution with that which would be 
tained were the concept of a single-component Hertz vector to be 
employed Now for the scalar problem the solution whose normal derivative 
vanishes on the sheet has been given by Macdonald (3) in a form equivalent 

H?)(kR cosh p) du — = H(kS cosh p) dp. 
Hs 
Let us suppose that this represents the Hertz vector for the solution of the 
ector problem. Identifying U with II,, we find 
E ‘| k?} H'‘- (k R cosh p) du 
i y~ 


‘| a ke] | HP (kScoshp) dp, (31) 
=\cy* P 
Hs 

yhich satisfies all the required boundary conditions except for that imposed 
t the edge of the sheet. Here (31) has a singularity of order r~! and it is 
this fact which compels its rejection as a possible solution of the physical 
ector) problem 

We now return to our solution (30). Whilst this is exact, for many pur- 
ses an approximation assuming /R, large compared with unity is suffi- 

Consider therefore the integral involving R. Since kR cosh pu kR, 

throughout the range of integration, we may replace the Hankel function 

the first term of its asy mptotic expansion to give 


yt 


: > : d 
HY R cosh p) dp ~ | ) im4 | e—tkR cosh B ‘ 
J J \7kR 4 vcosh 
LR 
Making the substitution 7 (2k R)sin $y, this becomes 
? , P e-' dr 
rf 1 ;, . — (32) 
N \n J v{(7?+kR)(7?+ 2k R)} 
k . 
nere Ty Zz | ‘0 Jeos 1(@ G,). 
\ R R, 7 
\n asymptotic expansion of (32) then follows by putting 7 equal to its lower 
it value in the non-exponential factor of the integrand. In this way we 
tain 
2 / 2 | 
| H‘ i R cosh p) dp ~~ ite a € ikR ind FY T ‘a (33) 


kl \7R\(R+R,)J a 
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va 


where F(u) ( e-ir? dy. 


Similarly 


Fea 2 /{ : \ .-ikS—inia BP 
(2) LS cos ~~ / , ___ \p-tkS—i1 f Te). ‘ 
J H\?(kS cosh p) du kal \7R (SERS (—Tg), (34 
where Ts, 2 I ia Joos 3(0+-6,). 
N \S- R, 


Hence 
2k 


E, ~ /( \ ikR,-—in Leos 1g cos 19 
' AN \7rR, 175 , z: 


oe? ke) /| - \, ikR+ini§ F( 7.) 
CYCY, kal \7R,( R+-R,)) 


o L /{ 2 er 

i. k2 / ( ikS+in + F( - ) (35 

; / 7 Tg)+ \v 
boa ka \7R,(S t R,)) 

This result, although approximate, ceases to be valid only for positions of 


both source and receiver within about a wave-length of the diffracting edge. 


5. Further remarks 

It is obvious from the method employed that our solution must satisf) 
the edge condition: for the field has been expressed as a double integral over 
plane waves each of which has the required behaviour at the diffracting edge. 
This may appear surprising in view of the double differentiation which our 
solution involves, and it is therefore worth while to derive the edge singu- 
larity. For convenience we use the approximate form (35) rather than (30 
this imposes little restriction upon the generality of our conclusions, demand- 


ing only that the source be more than a wave-length from the edge. A further 


simplification is to assume that z = 2», giving R, = r+-ry. For small rw 
then have 
R ~ 7—r cos(@—84,), S ~ ~—rcos(6+6,), 
TR & \(2kr)cos 4(6—8,), Tg ™ 4/(2kr)cos $(6+-4,). 


Carrying out the differentiations contained in equation (35) and retaining 
only the highest order terms in 1/r, 


: /? : e—ikr p—tkry { 2 ) 
BE, ~ k? | = etna cos 10) ( cos 164) cos 4, + (2—3 cos 8@,) 
A 7 \ (kr) kr, ales 2kry 


(36 


as r > 0. The edge singularity is thus of order r-?, which is in accordance 


with the accepted behaviour of the physical field. 
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One interesting feature of this result is the way in which the form of the 
ingularity depends upon the direction in which we approach the edge. In 
particular, if r > 0 in the direction 6 = 7, E, remains finite there. An 
explanation, however, is not far toseek. For the singularity may be regarded 
1s that due to a line source coincident with the edge of the sheet and of 

wiable angular spectrum, @ 7 being a direction of zero intensity. It is 
this secondary source which is responsible for the luminous nature of the 
liffracting edge. 

Moreover, this behaviour is not restricted to the diffraction of a spherical 
wave. If the component of the electric field normal to the sheet is calculated 
for the diffraction of a plane or cylindrical wave, the same dependence on 6 
s found, the singularity vanishing if r > 0 in the direction @ = z. But in 
these latter cases the behaviour is less conspicuous since it is customary to 
work in terms of that component of either E or H which is parallel to the 
edge 

Finally it may be remarked that the treatment adopted in this paper is 
equally applicable to the determination of field components other than E,,. 
Thus, for the complete field we find 


E a" 14 cos 40, H2(kR,) +- = ee on 
(7% 2 Cx0y, 2 CXCYy 
E ee cos 44 cos Ld, HY (kR,) + ‘| ai i \Ip— ( ae + kN, 
(77) 7 ‘ 2\CyoYy, *  2leyeyy 
EB i Cc I, ” o* 1. 
2 oz Yo 2 < CYUg 
kY (z—z . KY of, kY ol, 
H . . 0) cos $4 cos 40, HY (kR,) a aR - fh 
RR ty) = CZ a Ox 
‘ kY | en) « 9 
H ©’ sin 36 cos 36, HP (kR,), 
Ry (779) : - 
H I el, ky ol. 
_ ari) 2 OXy 
vhere 
[, ( HP (kReoshu)du and I, [ H‘?(kS cosh p) du. 
Hs 


T 
| 


the presence of a non-zero component H, should be particularly noticed. 
The author wishes to express his thanks to Dr. P. C. Clemmow for many 
helpful suggestions; also, to acknowledge receipt of a D.S.I.R. grant for the 


period during which this work was completed. 
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THE EXTERNAL MAGNETIC FIELD OF A SINGLE THICK 
SEMI-INFINITE PARALLEL PLATE TERMINATED BY A 
CONVEX SEMI-CIRCULAR CYLINDER 
By N. DAVY and N. H. LANGTON (The University, Nottingham) 
[Received 25 September 1951] 


SUMMARY 


fhe theory of the external field of a common type of electromagnet polepiece is 


lerived. This is a single plane slab of magnetized ferromagnetic material with a 
nvex, exactly semicircular cylindrical end, i.e. shaped like a book-cover with a con- 
ex spine. Num«e rical tables and a graph are given to illustrate the variation of the 
field strength along principal lines. Results which may be of interest in molecular 


eam work are included. 


1, Introduction 

\ comMON type of electromagnet polepiece consists of a plane slab of mild 
steel with a convex rounded end, shaped like a book-cover with a semi- 
ircular convex spine (see Fig. 1). It is of interest to derive the exact theory 
of the field of a complete magnet having two such opposite polepieces, with 
in air gap between them. Cockcroft (1) gave an approximate theory. 
Though the exact theory for two polepieces seems intractable at present, 
the exact theory of the field of one such polepiece is attempted in this paper. 
Page (2) solved the case when the end profile curve was a cycloid. The 
electric field of a single charged conductor of a similar shape may also be 
finterest in connexion with the mechanism of lightning discharges. The 
iathematical calculations for two-dimensional systems are the same in 
both the magnetic and electric cases. It is proposed, then, to determine the 
two-dimensional external magnetic field of a single semi-infinite body of the 
shape shown in Fig. 1, that is, of a polepiece BGODB’ whose shape is that 
fa plate, with parallel sides extending to infinity in one direction and 
terminated by a convex semicircular cylinder in the other. The edge 
BGODB’ is regarded as a magnetic equipotential surface. This assumption 
8 justified by its successful application to other magnets, and entails the 
ondition that the material used shall be ferromagnetic and in the un- 
saturated state. Although BGODB’ consists only of two straight lines and 
i semicircle, its field cannot be obtained by any of the three well-known 
ransformations W A/z, Wi Alnz,and W = A In(z—a)/(z+a) involv- 
ng circles, or by Richmond’s method. A different method of attack is 
hecessary. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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IZ 
Fic. 1. The z-plane. Cross- 
section of a polepiece or 
conductor. 


2. The transformations 

(a) The geometrical transformation. The edge BGODB’ in the z-plane can 
be regarded as a curvilinear triangle bounding the external region, so that 
the angles at the corners G, D, and B are z, 7, and 27 radians, respectively. 
It can be transformed into the real axis of another diagram, the c-plane 
(not shown), and its exterior region into the upper half of that plane bya 
method originally due to Schwarz (3) which has also been used by Daymond 
and Hodgkinson (4), in a certain aerofoil problem. The geometrical trans- 
formation gives z as a function of two solutions of the hypergeometri 
equation, the independent variable being c. Here ¢ is both the variable of 
the new plane and the square of the modulus k of the elliptic integrals f, 
kK’, E, and E’ in terms of which z is expressed. At D let k? = c¢ = +1 
at O let c +2: at Gletc +-oo; and at B and B’ let c = 0. 

Quoting from Daymond and Hodgkinson (4), the appropriate geometrica 
transformation, expressed in the present notation, is 

Z ia{(2—c) BE’ —cK'}/{(2—c) E—2(1—c)K}+-a(1+2)/2. (I 

Here the origin of coordinates is the point O, the ‘nose’ of the body, as 
the full width, i.e. twice the radius of the end, and 7 /(—1). 

(b) The magnetic or electric transformation (see Fig. 2). It is now necessary 
to effect a second transformation, the magnetic transformation, to link the 
magnetic properties of the field of the body in Fig. 1 to its geometrit 
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Let the potential of the body BGODB’ be 
This body extends to infinity at B and B’, where the variable c is equal to 
ero. The potential V of all points outside BGODB’ in Fig. 1 is less than 

J, and may be supposed to range downwards to minus infinity in value. 


properties +-V, everywhere. 





onsider a diagram (Fig. 2) called the W-plane, where W = U+iV. Here 
b G O D B 
Fic. 2. The W-plane 
nfinitely long horizontal straight line BGODB’ has the equation W U+iy,. 
represents the magnetic state of the polepiece in Fig. 1, which is at constant 
ial J Points below BGODB’ correspond to points external to the polepiece 


in Fig. 1. 


isthe magnetic potential already mentioned and U is the stream function, 
mstant along a line of force. 
vy = U-+iv,, extending to infinity at both ends, represents the values of W 
mo the edge of the body BGODB’ in Fig. l. 
ely real values as along BGODB’, withe 


The region in Fig. 2 extending below the horizontal line W = U-+ VY, to an 


In this figure the horizontal straight line 


Along it c has the same 
0 at the two ends at infinity. 


nfinite depth corresponds in its values of U and V to points external to 
BGODB’. We wish to transform the line W 


xis of the c-plane, and the region below it extending to infinity into the 


U--<.V, in Fig. 2 into the real 
0 > 


ipper half of the c-plane, point for point. This is done by first writing 
W—wW, W, 


4 a new variable, and then letting W, = A/c where A is a con- 


tant, so that finally 
W —w, 


nd dv dc 


A c | 
Alc? } 
rhe first process, in effect, transforms the whole region extending below the 


ne W iV, 


) in Fig. 2 into the lower half of a W,-plane, below the real 
xis of that plane, and the second process is a complex inversion of that 
wer half of the W,-plane about a point in the real axis, transforming it into 

the upper half of the c-plane. The equation 

W = A/e+uM (3) 

iorms the magnetic transformation. Combining equations (3) and (1) so as 
eliminate c we express W as a function of z. The magnetic pole strength 
er unit depth residing along any length of the edge BGODB’ in Fig. 1 is 
siven by the modulus of (W,—W,)/47, where W,, W, are the values of W at 
the ends of the stretch considered. Thus, between the points D and G in 


Nig. 1, points at which c | and 00, respectively, the pole strength per 
init depth is A /47. The physical meaning of A is 47 times the pole strength 


per unit depth on the rounded end of the body in Fig. 1. 








































118 N. DAVY AND N. H. LANGTON 
3. The magnetic intensity 


The magnetic intensity (field strength) F at any point in the 2-plane 
aiencdicsias p_|aW| _ \aw faz 
7 | dz | | de} de} 


The external medium is assumed to be a vacuum. By (2) dW /dc is equal 


(4 








to —A/c*; and dz/dc is obtained by differentiating the value of z given by 
(1). This is done by the aid of differentials of K, K’, E, and E’ (Cayley’s 
Elliptic Functions, p. 48), and we find 


dz/de = |—3tamc/4{(2—c) EH —2c’ K}?|. (5 
Hence F = |dW/dz| = |41A{(2—c)E—2c’K}?/3azc}|. (6 


We now need the numerical values of c, K, K’, EZ, and E’ at the points 
where F is to be calculated. Points along the edge BGO or B’ DO and along 
the axis OZ are of interest in this connexion. At the point D, where c = 1, 
F - | - 41A/38an| = 4A /3a7 (7 
since A is real. 
(a) Calculation of F and z along B'D 
In this case 0 < ¢ < +1, and there are no special difficulties. Values of 
K, K', E, and E’, obtained from tables of elliptic integrals, are substituted 
in equations (1) and (7). At the end points of this length, viz. at B’, we have 
c = 0, and at D,c = +1. 
(b) Calculation of F and z along the edge of the quadrant DO 
Here +1 <c < +2 and no corresponding tables of elliptic integrals 
seem to be available. We have to use the fact that when k becomes 1/k and 
c becomes l|/c, then K becomes k(K-+iKk’'), K’ becomes kK’, E becomes 
{E—c'K—i(E’—ck')\/k, and E’ becomes E’'/k. Also, writing i/¢ = C, 
along DO we get 
z = ia{(20—1) EB’ —CK"\/{(2C—1) E+ C’K —i(2C—1) EB’ +iCK}+ 
+a(1+2)/2, (8) 
F = |—4iA{((2C—1)E+C’K—i(2C—1)E’+iCK"}2/3az|. (9 
In this case K, K’, E, and E’ are now functions of C and also +} < C < +1. 
Hence their values are obtained from the usual tables. 
(c) Calculation of F and z along OZ 
Along OZ, except at O itself, c must have a complex value, since its real 
values correspond to points along the edge BDOGB’. Let 


¢= 2/(l+c,) or c, = 2/c—1. (10) 


3y this transformation, we have introduced a variable c, which has wholly 
+o, along the axis OZ. 


imaginary values, that isc. = —ia where 0 < a < 
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At the nose O where c 2, we have c, = 0. At the point D, where 

1, we have cy, 1. At the point G, where c = 00, c, = —1, and 
it the point B, where c 0, we have c, = oo. At the point Z, where c = 0, 
we have Cy ico. The introduction of c, has also introduced symmetry 


about the point O. Also 
c 2/(1+-¢,) 2/(1—ia) = 2/(1+a?)+72a/(1+.0?). 


Writing x = sinA/(1+-cosA) = tan(A/2) 
we get c 1+cosA+isinA = 1+e!, (11) 
Co i tan(A/2). (12) 
Although tables of K, K’, Z, and E’ for values of c, = —ia do not seem 


to be availak'e, certain tables published by Cambi (5) give some values 
of K and K’ (but not of EZ’ and £) when c = 1+ e!4. However, E can be 
calculated by a method based on the use of a special function tabulated by 
Cambi (Table VI, p. 244, loc. cit.). 
Finally, values of E’ are calculated by Legendre’s equation, 
KK'+ E'K—KK’ = a/2, 

which still holds even when c = k? is complex. Numerical values of z and F 
are obtained by substituting the values of K, K’, EF, and £’ in equations 
1) and (6). These laborious calculations have to be carried out afresh for 
each point on the line OZ. The labour would be much reduced if tables of 
withmetical values of K, K’, E, and E’ as complex functions of the complex 


> 


variable k? = c = p+iq were available. 


4. Results 

The values of the field strengths at various points along B’D, DO, and 
0Z are given by Table I and represented graphically in Fig. 3. Ordinates 
are proportional to field strengths, abscissae to distances measured along 
the three sections of the line B’DOZ in Fig. 1. As expected, there is a 
maximum of field strength and surface density at the ‘nose’ O, contrasting 
with Page’s eycloidal tipped polepiece for which F and o are constant over 
the end. There is a point of inflexion in the graph at the corner D, corre- 
sponding to the geometrical ‘kink’ there. A further result, not revealed by 
Fig. 3, is that along the line OZ the values of log F and the variable A, where 

|+-e', are connected by a linear relation over the initial range starting 

Irom Q, 

|. Field strength F on the axis OZ (see Figs. 1 and 3). The field strength 
Fe.m.u. is expressed as a multiple of the constant 4.4 /3az of equation (6), 
&, Table I gives 3a7F'/4A. The distance z of the field point from O is 
expressed as a multiple of —ia/2, i.e. the table gives —2z/ia. The modulus c 
is complex, except at the point O, where c = + 2-000. 
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TABLE | (see also Fig. 3) 


{ I°Q51I { 18090 { 1°5878 {| I*3090 
é 2000 ; 4 : 

(+2:3090 =|+7°5878 |+72-8090 = | +2-9511 
— 2z/ta 0:000 0°254 0°556 0°982 1°590 
3amF/4A 1°719 1°394 1*134 0°930 0°753 





| 1'000 
{+72 1-000 
2°452 
0600 


2. Field strength F e.m.u. on the circular edge OD (see Figs. 1 and 3). Here 


s is the length of arc measured along the edge from O t« 
radius is a/2 and Table IT gives 2s 
At the point D, s 


a, i.e. length of are as 
a/2xa/2 = 1-571a/2. The modulus c 


than unity, except at D, where it is equal to unity, and F 


02 
ToZ | The Nose 

1 

12 O08 O04 O 4 

Fic. 3. 

line ZODB’ in Fig. 1. 

distances, in multiples of the radius a/2, 

measured from the nose (1) to the right 

along the edge ODB’, (2) to the left, along 
the axis of symmetry OZ in Fig. 1. 


he Point 


Field strengths at points along the 


The abscissae are 


TABLE II 


2°000 1°704 1°490 1°333 1218 1°133 1'072 
2s/a 0000 «0239 0°475 0°703 o'919 I'lig 1°296 
3am F/4A 1°719 1°705 1°662 I°591 1°500 1°387 1°272 


) any point. The 
a multiple of a/2. 





is real and greater 
4A /3ar. 





lo B’ 


— 


1°03I 1°005 1°000 
1440 1°530 1°57! 
1°156 1°057 1000 


3. Field strength F e.m.u on the straight side DB’ (see Figs. 1 and 3). Here 
s is the distance of the point concerned from D along the straight edge. 


Table IIT gives 2s/a, i.e. distances as a multiple of a/2. 
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The modulus ¢ is real and less than unity except at D. 


TABLE III 


0°933 0553 o'd21 0°750 0"500 0250 
0°404 os160 1*505 2°O15 12°79 50°00 
7 0°745 0°037 0°542 0°450 0°230 0°0999 


Two results likely to be useful in magnetic and molecular beam measure- 


ments, when BGODB' isa polepiece, are obtained by differentiating (6), viz. : 


dF /dz 16.4$(2—c)K —2E¥(2—c) E—2c’K}3/3a?xc'|, (13) 
dF2/d2 = |128iA%{(2—c)K —2E}{(2—c) E—2c’ K}5/9a3nc8|. (14) 


The exact internal field of a charged conductor of the same shape as 
BGODB’ has already been calculated (6). 
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SOLUTION OF ORDINARY LINEAR DIFFERENTIAL 
EQUATIONS WITH VARIABLE COEFFICIENTS BY 
IMPULSIVE ADMITTANCES 


By W. J. DUNCAN 


(Department of Aeronautics and Fluid Mechanics, the University, Glasgow) 
[Received 15 May 1952] 


SUMMARY 
Linear differential equations of any order and sets of such equations can be solved 
by means of impulsive admittances and the solution expressed in terms of a simple 
definite integral or a sum of such integrals. Hitherto the method has been applied 
only to equations with constant coefficients, but the method is now extended to 
equations and sets with variable coefficients. The process of solution, suggested by 
a physical analogy, is simple and direct. 


1. Introduction 

THE method of impulsive admittances is suggested by a physical problem. 
Suppose that the differential equation considered is linear and represents 
the motion of a body or a system under the action of a force which varies 
with time in a prescribed manner. Then the differential equation consists 
of a dynamical variable operated upon by a linear differential operator and 
equated to the known function of time which represents the applied force. 
This force may be supposed to be built up from instantaneous impulses, and 
if we know the free motion which follows such an impulse we can obtain 
the required solution as a definite integral. The function which represents 
the free motion of the system, supposed to be initially at rest in its datum 
state and following a unit impulse, is defined as the impulsive admittance. 
This is a species of Green’s function. 

Cauchy developed a method of solving linear ordinary differential 
equations with constant or variable coefficients which is essentially the 
same as that of impulsive admittances and an example of his method is 
given by Temple and Bickley (1); reference may also be made to Goursat (2). 
However, Cauchy’s method does not appear to be very well known and it 
has received comparatively little technical application. It seems desirable 
therefore to give a new exposition in a very general form and with emphasis 
on the simple physical interpretation of the method. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 1 (1953)] 
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The method of impulsive admittances has hitherto been explicitly 
developed for linear differential equations of the first and second orders 
with constant coefficients and for sets of such equations in several variables 
3, 4, 5); such sets may be supposed to be the dynamical equations of a 
system with several degrees of freedom. In this case the impulsive 
.dmittances are functions of the single variable (t—r), where t¢ is the time 
it which the dynamical variable is measured while 7 is the earlier time 
it which the unit impulse was applied. When the coefficients in the 
differential equations are functions of the independent variable (time) the 
impulsive admittances clearly become functions of the two variables ¢ and 7, 
for the response to a unit impulse now depends on the time at which it is 
pplied. In the present paper this method is applied systematically to linear 
differential equations of any order with variable coefficients and to sets of 
such equations. tT 

Any solution obtainable by the method of impulsive admittances could 
so be derived by the method of the variation of parameters (7, 8). How- 
ever, the method of impulsive admittances is far more direct and simple in 
ipplication. The method of the variation of parameters may fairly be 
described as an ingenious mathematical trick which happens to be success- 
ful, whereas the method of impulsive admittances is based on a simple 
physical idea which suggests in detail the expression for the solution. 

It is assumed throughout that the conditions for the validity of the 
differentiation of the definite integrals which appear in the investigation 
re satisfied 

The writer is well aware that the theory given in sections 3 and 4 below 
could be stated more concisely by means of the matrix notation, but he has 
not adopted this method of presentation in order to avoid giving the 
impression that the method of impulsive admittances is specially bound 


up with matrices 


2. Single equation of the nth order 
The equation to be solved is 


ad" Jn 1 
P (t) u >» ey! 


qd , 
? dt" | n—-1' T 


dt" 1 vet Pi(t)q Q(t), (2.1) 


where P(t) and Q(t) are given functions of the independent variable ¢. The 
impulsive admittance a(¢,7) is a function of the two variables t, 7 and is 
defined as follows: 

\n approximate method for carrying out response calculations has been developed by 


Tustin (6) based on the use of a ‘Delta input’, i.e. a force whose graph is an isosceles triangle. 
In the limit when the base of the triangle tends to zero the response to this becomes effec- 


tively identical with an impulsive admittance. The admittance method has also been 


} 


eveloped in America under the name of the ‘pulse method’ (7). 
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(a) When 7 is constant and t > 7, a(t, 7) is a solution of the homogeneous 
equation 
dq d” lg : 

P (t)—=+- P. _.(#) = +....4- Bit)a = 0. 2.) 

n ) an l n 1 ) dt” a. 1 of \q ( 4 
(b) usr 7) = © 
x (7,7) 0 

(2.3) 


a(t, r 
where v(t,7) =< a ) (2.4) 
cot’ ; 
The function a(f,7) is not defined for ¢ < 7 and the differential coefficients 


which appear in (2.3) are defined as limits for t approaching 7 from above. 
A solution of equation (2.1) is then given by 


t 


J) = [ Q(r) a(t, 7) dr. (2.5) 


Moreover, f(t) satisfies the boundary conditions 
f(0) 0 
f'(0) 0 


, (2.6) 
f"-(0) 0 | 
where f"(t) = d’f/dt’. 
The proof is the same for equations of all orders and it will suffice to take 
n = 3. We have from (2.5), by the rules for the differentiation of definite 
integrals, 


t t 
f'() = Q(t) aft, t) 4 [ Q(r) a’ (t, 7) dr [ Q(z) x’ (t, 7) dr (2.7) 
0 0 
t 
by (2.3). Similarly {"@) = O(r) x" (t, 7) dr (2.8) 
0 


of, 
while 
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Hence 
Pt) f"(HW+PR(OS"(OW+ROSf' O+ROfO 


Q(t)Pa(t) «"(t,t)+ | Q(r) F(t,7) dr, (2.10) 
0 
where 
F(t, 7) P,(t) x” (t, 7) +- P(t) x” (t, 7) + P,(t) x’ (t, 7) + Po(t) a(t, 7) = 0 
by the defining property (a) of the impulsive admittance. Then, on account 
of the last of equations (2.3), the expression on the right-hand side of 
2.10) reduces to Q(t). Thus f(¢) satisfies (2.1) and by equations (2.5), (2.7), 


ind (2.8) 


f(0) = f*(0) = f"(0) = 0 
in aceordance with (2.6). The complete primitive is obtained by adding the 


complementary function to f(t). 


3. Set of m equations of the first order 


The set of equations to be solved is typified by 


S444 B..(t)q.) = Q,(t). (3.1) 
=A dt 7 J 


Then the typical impulsive admittance «,,(t,7) has the following defining 
| kl g 2 
properties 
a) q x,(t,7) (2 constant) 


is a solution of the set of homogeneous equations obtained from (3.1) by 


making the functions Q(t) vanish. 


b) > A.,.(7) a7, T) on 
1 
where bn | 
ind 0,7 0 (r +f, 
while the determinant |A..(7)| must not vanish. 


The complete solution is then given by 


t 
q(t) > Q (7) Xy,(t, 7) dr. (3.2) 


p=1 4 


This can be proved as follows. Equation (3.2) yields 


la (t) m é 
J SQ (t) ve) (t,t) 4 > | Q@ (7) Xgp(t, 7) dr. 
1: 


y = t 
dt Ses . 


0 
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Therefore 


om 7 m. 
S14, (t) 4 B(tg,| = 3 t)A,.(t 
a > vil dt ra\ Ms) = p 3 Qy( t). rs( )a Xsp| »t) 


1 p s 
[ Q, (7) F(t,r) dr, (3.3 
siti 
where F(t,7) Ss {A,.(t) x Xen (ts 7) + B,,(t) Xsn(t, 7)} @ 


s=1 

by the defining property (a) of the impulsive admittances. Also, by the 
property (5), the first summation on the right-hand side of equation (3.3) 
reduces to Q,(t). Hence the set (3.2) satisfies equation (3.1). On account 
of (3.2) the boundary conditions are 


q(0) = 0. (3.4) 


4. Set of m equations of the second or higher order 
The set of second-order equations to be solved is typified by 


m 


> vl () 24+ B, st) B+ Catto. — Q,(t). (4.1) 


The impulsive admittances now have the defining properties: 
(a) Uz = %,(t,7) (2 constant) 
is a solution of the set of homogeneous equations which is a special case 
of (4.1), 
(b) Xp4(T, T) 0, 
m P 
(c) 7 A,,(7) og(7,7) = 5,, 
s=] 
where 5,,is Kronecker’s symbol as before, and the determinant | A,,(7)| does 
not vanish. The solution is again expressed by equation (3.2). 
We now derive by use of the property (b) of the admittances that 


t 
_ he | Q,,(7) a(t, 7) dr, (4.2) 
while Ss Op (t, t) >) f Q (7) xg,(t, 7) dr (4.3) 
dt? p= “ep p= 1 5 J sa . 
Accordingly 


m 


2 | A,. g@ “Os LB. (t) +C,,(t)q,! : Y ys Q,(t)A, 4 (t) v, p(t, t)-+ 
ie dt | p=1s=1 


. 


5 | Q(z) F(t, 7) dr, 
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where 


rs(l) Ys (0,7) 1 B,,(t) Xgp(t, 7) tC, (t) Xsy(t,7)} = 


rs 


Fit.r) > {A 


by the defining property (a) of the admittances. Also, by the property (c), 

the first summation on the right-hand side of equation (4.4) reduces to 

Q,(t); this verifies the solution. From (3.2) and (4.2) it follows that the 

boundary conditions are 

q.(0) = 0 r 

: |, (4.5) 

q7.(0) 0 | 

The foregoing method can obviously be extended to a set of m equations 

f the nth order and the solution is still expressed by (3.2). The admittances 
have the property (a) as above, while condition (5) is amplified to 


¥(7,7) = 0 


see equation (2.4)). 
Condition (c) becomes 
> A,,(7) a (7,7) = 3y, 
s=1 
where the functions A,.(t) are the coefficients of the differential coefficients 
of order n and A..(7) does not vanish. 
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HALVING THE INTERVAL IN A TABLE WHEN 
FIRST DERIVATIVES ARE GIVEN 


By HAROLD JEFFREYS (St. John’s College, Cambridge) 


[Received 18 October 1951] 


I RECENTLY had to evaluate numerically a number of integrals of the form 
1 
| x" F(x)G(x) dx, 
0 
where F and G@ were given at rather wide intervals, so that the Gregory 
formula would be very inaccurate. On the other hand, F’(x) and @'(a) 
were given for all tabular values of x, including the end points. In these 
conditions the Euler—Maclaurin formula, including first derivatives, would 
be applicable, but the actual evaluation of the derivative of the integrand 
would be rather troublesome. 
The method actually adopted was to interpolate F(x) and G(x) to the 
mid-points of the tabular intervals by the formula 
f (3h) = HS (O)+AS(A)}— Shih (A) '()}- (1) 
It is easy to verify that this is exact for a cubic polynomial. For f(x) =# 
at interval 1, it gives f(}) = 0. Bessel’s formula, including second differ- 
ences, would give — 4; the error of (1) is therefore a ninth of that of the very 
accurate. Bessel formula even when all the information needed to apply 
the latter is available. 
When the interval was halved in this way the third differences were 
small enough for the Gregory integration formula to be satisfactory. 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 1 (1953)] 
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